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To  date,  very  little  research  on  the  dynamics  of  parallel  kinematic 
mechanisms/machines  (PKMs)  exists.  Of  the  work  available  in  the  literature,  most  have 
been  for  generating  the  equations  of  motion,  inverse  dynamics,  and  structural  dynamics. 
Insight  to  the  structure  and  properties  of  PKM  dynamics,  effects  of  strut  inertias, 
dynamics  interactions  with  various  kinematics  error  sources,  formulation  of  dynamic 
performance  indices,  and  dynamic  balancing  for  PKMs  have  only  recently  been  provided. 
Due  to  the  above  research  or  lack  thereof  has  created  an  increasing  awareness  of  the 
dynamic  needs.  Therefore,  the  thrust  of  the  work  presented  in  this  dissertation  is  to 
establish  a dynamic  parameter  identification  approach  for  parallel  kinematic  machines. 
This  approach  is  to  acquire  parameter  values  necessary  for  modeling  the  dynamics  of 
parallel  kinematic  machines,  e.g.,  mass,  inertia,  and  friction. 

To  determine  the  best  approach,  several  identification  methods  were  compared  and 
evaluated  in  terms  of  accuracy,  convergence,  ability  to  handle  noisy  data,  and  ease  of 
implementation.  The  main  comparison  studies  were  conducted  between  the  least  squares 


xv 


(LS)  method  and  the  unscented  Kalman  filter  (UKF)  method.  To  illustrate  the  superiority 
of  the  unscented  Kalman  filter  (method  of  preference),  both  the  LS  and  UKF  methods 
were  applied  to  the  inverted  double  pendulum  case.  Using  results  found  in  the  literature, 
the  extended  Kalman  filter  capabilities  were  compared  to  these  two  methods,  further 
substantiating  the  superiority  of  UKF  method  for  nonlinear  systems.  The  comparisons 
were  simulation-based.  In  addition,  force-based  and  energy-based  modeling  methods  are 
compared  to  determine  if  there  is  any  benefit  other  than  reduced  modeling  effort  in 
deriving  the  theoretical  model.  A significant  improvement  using  energy-based  method 
was  not  generally  realized.  To  experimentally  validate  the  capability  and  the 
implementation  of  the  UKF  method  for  dynamic  parameter  identification,  the  University 
of  Florida  Space,  Automation  and  Manufacturing  Mechanisms  Laboratory’s  PKM 
dynamic  parameters  were  identified. 

The  simulated  and  experimentally  validated  results  show  that  the  unscented 
Kalman  filter  performs  well  in  identifying  the  system  parameters,  eliminating  the  need 
for  linearization.  This  merit  corresponds  to  the  significant  reduction  in  modeling  effort 
required  to  generate  the  basic  algorithm  used  for  the  system  identification.  The  overall 
impact  to  technology  is  provision  of  a parameter  identification  method  that  yields 
improved  rigid  body  dynamics  models  of  parallel  kinematic  mechanisms.  This  will  lead 
to  the  ability  to  introduce  advanced  controllers  that  exploit  the  model  dynamics  to 
improve  system  performance.  In  addition,  the  presented  method  sets  the  foundation  for 
formulating  an  online  tuning  algorithm. 
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NOMENCLATURE 


(a„)_  absolute  acceleration  of  the  center  of  mass  of  the  body  ‘B’  of  strut  in  inertial  frame 
{W}  coordinates. 

a t absolute  acceleration  of  the  center  of  the  mass  of  the  moving  platform  in  inertial 
frame  {W}  coordinates. 

(afi).  absolute  acceleration  of  the  center  of  mass  of  the  body  ‘E’  of  strut  in  inertial  frame 
{W}  coordinates. 

bj  base  platform-strut  joint  for  joint  i;  i=l  ,6. 

d k measured  (desired)  output 
det[  ] determinant  of  matrix  [ ] 

reference  residual  i. 

F double  pendulum  example’s  controller  input  force  (Chapter  3) 

F„  6-component  vector  [fU|  ,F(h  ,---Faf  ]'  of  the  six  actuator  forces ) , i=l,6  [units: 


(N,  N,  N,  N,  N,  N)t] 

{F}  coordinate  system  (XF,  YF,  ZF):  frames  with  the  same  orientation  as  frame  {W}  but 


with  origin  located  at  the  center  of  mass  of  the  moving  platform 


F 

c 

T 


cen 


cen 


to  J Ffl  wrench  due  to  the  centripetal  type  accelerations  of  the  strut  masses  [units: 
(N,  N,  N,  N-m,  N-m,  N-m)T],  defined  in  the  inertial  frame  {W}. 


6-component  vector  of  the  (forces,  moments)  contribution 


F 

cor 

T 

_ cor 


to  J F„  wrench  due  to  the  Coriolis  type  accelerations  of  the  strut  masses  [units:  (N, 
N,  N,  N-m,  N-m,  N-m)T],  defined  in  the  inertial  frame  {W}. 


1 

£ 

i 

- 

F r 

T rl 

i 

1 

ex  t 

ext  J 

6-component  vector  of  the  (forces,  moments)  due  to  the 


externally  applied  forces  and  moments  acting  on  the  PKM  [units:  (N,  N,  N,  N-m, 
N-m,  N-m)T],  defined  in  the  inertial  frame  {W}. 


6-component  vector  of  the  (forces,  moments)  contribution 


to  J F„  wrench  due  to  the  joint  friction  [units:  (N,  N,  N,  N-m,  N-m,  N-m)T], 
defined  in  the  inertial  frame  {W}. 


6-component  vector  of  the  (forces,  moments)  contribution 


to  J F„  wrench  due  to  gravitational  acceleration  acting  on  the  platform  mass  [units: 
(N,  N,  N,  N-m,  N-m,  N-m)T],  defined  in  the  inertial  frame  {W}. 


6-component  vector  of  the  (forces,  moments)  contribution 


to  J Ffl  wrench  due  to  gravitational  acceleration  acting  on  the  strut  masses  [units: 
(N,  N,  N,  N-m,  N-m,  N-m)T],  defined  in  the  inertial  frame  {W}. 
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[f  r T T] 

L pm  pm  J 


6-component  vector  of  the  (forces,  moments)  contribution 


pm 


pm 


to  J F„  wrench  due  to  the  linear  and  angular  accelerations  of  the  platform  mass 
[units:  (N,  N,  N,  N-m,  N-m,  N-m)1],  defined  in  the  inertial  frame  {W}. 


6-component  vector  of  the  (forces,  moments)  contribution 


to  J Ffl  wrench  due  to  the  linear  and  angular  accelerations  of  the  strut  mass  [units: 
(N,  N,  N,  N-m,  N-m,  N-m)T],  defined  in  the  inertial  frame  {W}. 


6-component  vector  of  the  (forces,  moments)  contribution 


to  J Fa  wrench  due  to  un-modeled  dynamics  [units:  (N,  N,  N,  N-m,  N-m,  N-m)T], 
defined  in  the  inertial  frame  {W}. 

(F  ) magnitude  and  sense  of  the  translational  friction  force  for  friction  in  strut  i, 


direction  along  (s3). ; i=l ,6. 


[I]  a 3x3  identity  matrix 

(Ixx),  (Iyy),  and  (IZ2)  principal  mass  moment  of  inertias  of  the  platform  (kg-m2),  defined 
about  the  center  of  mass  of  platform  with  principal  axes  |s2,s3,s23  I rrespectively. 

(Ixxb)i,  (Iyyb)i,  and  (IZZb)i  principal  mass  moment  of  inertias  of  the  nonextendable 

portion  of  the  strut  (kg-m2),  defined  about  its  center  of  mass  with  principal  axes 
{s 2, s 3, s23 Irrespectively.  One  set  for  each  strut  i. 

(Note,  lyyb  is  moment  of  inertia  about  s3  axis  and  hxb^zzb  are  about  the 
corresponding  orthogonal  axis  to  § ) 
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(Ixxe)i,  (lyye)i,  and  (Izze)i  principal  mass  moment  of  inertias  of  the  extendable  portion 

of  the  strut  (kg-m2)  defined  about  its  center  of  mass  with  principal  axes 
{s 2 , s3 , s 23 }.  respectively.  One  set  for  each  strut  i. 

(Note,  Iyye  is  moment  of  inertia  about  s3  axis  and  IXXe=Izze  are  about  the 
corresponding  orthogonal  axis  to  s3 ) 

J,  J"1  (6x6)  pseudo-Jacobian  matrix  and  inverse  of  the  pseudo-Jacobian  matrix, 
respectively 

JLS  k Jacobian  matrix  used  in  the  least  squares  method 

(Zs3)i  vector  denoting  ith  strut  length  i along  strut  axis  (s3)i;  one  for  each  strut;  i=l,6. 

( Zs3  )i  vector  denoting  the  time  rate  of  change  of  the  strut  length  (l  )i  along  strut  axis 
(s3)i;  one  for  each  strut;  i=l,6. 

(Zs3  )i  acceleration  component  of  strut  that  is  aligned  with  unit  vector  (s3)j;  one  for  each 
strut;  i=l,6,  a vector  denoting  the  2nd  derivative  of  i 

(Lb s3).  a 3 -component  vector  denoting  the  location  of  the  center  of  mass  of  the 
composite  nonextendable  part  of  strut  i from  center  of  the  base  platform-strut  joint 
bj  along  (s3).  axis;  i=l ,6. 

(Lb  s3 ).  time  rate  of  change  of  the  distance  (l  ) along  strut  axis  (s3 ). ; i=l  ,6. 

Note,  Lb  =0  since  the  Lu  does  not  change  in  time. 

[Lb\  time  rate  of  change  of  LB , Note,  LB  =0  since  the  LB  does  not  change  in  time 
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(Z£s3 ),.  a 3 -component  vector  denoting  the  location  of  the  center  of  mass  of  the 

composite  extendable  part  of  strut  i measured  from  the  center  of  the  base  joint  bj 
along  (s3).  axis;  i=l,6. 

(Le  s , ).  time  rate  of  change  of  (i  f ) along  strut  axis  (s3 ). ; i=l  ,6. 

(l  f ) time  rate  of  change  of  (i  f ) ; i=l,6. 

/,  ,w,  length  and  mass  of  the  bottom  pendulum  (link  1),  respectively  (Chapter  3) 

1 2 , m2  length  and  mass  of  the  top  pendulum  (link  2),  respectively  (Chapter  3) 

M mass  of  the  cart  (Chapter  3) 

[M  p ] mass  matrix  relating  to  [f  J ,T  ' ] 

(mp)j,  (Ip)j  lumped  mass  and  inertia  tensor  of  the  moving  platform,  respectively. 
Center  of  mass  located  at  origin  of  {P} frame 

(mb)i,  (Ib)i  lumped  mass  and  inertia  tensor  of  the  nonextendable  portion  of  the  strut  i; 
i=l,6,  respectively. 

(me)i,  (Ie)i  lumped  mass  and  inertia  tensor  of  the  extendable  portion  of  the  strut  i; 
i=l,6,  respectively. 

{P}  coordinate  system  (Xp,  Yp,  Zp):  frame  rigidly  attached  to  the  moving  platform  and 
used  to  define  the  absolute  orientation  of  the  moving  platform,  origin  is  located  of 
the  platform’s  center  of  mass. 

Pi  moving  platform-strut  joint  for  strut  i;  i=l  ,6. 

Q 6-compenent  vector  of  the  generalized  coordinates  of  the 

system,  the  first  3-components  are  the  {W}  coordinates  of  the  absolute  position  of 
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the  moving  platform’s  center  of  mass  and  the  last  3 components  are  the  Euler 
angles  defining  the  absolute  orientation. 

q = \xc,  y c,z  c,<j>  ,0  ,y/~\  6-component  vector  of  the  generalized  velocity,  first 

derivative  with  respect  to  time  of  the  generalized  coordinates 
w r = (x  , yc  ,zc)7  3-component  position  vector  of  the  center  of  mass  of  the  platform 

specified  in  {W}  coordinate  system  [units:  (m,  m,  m)T]. 

(w  rt  )i  3 -component  position  vector  of  the  ith  base  -strut  joint  location  (bj)  defined 
relative  to  the  {W}  origin;  i=l ,6. 

("  rB  )i  3-component  position  vector  of  the  center  of  mass  of  the  nonextendable  part  of 
the  ith  strut  defined  relative  to  the  {W}  origin;  i=l,6. 

( " r )i  3 -component  position  vector  of  the  ith  moving  platform-strut  joint  location  (Pi) 

defined  relative  to  the  {W}  origin;  i=l,6 

("  rE  )i  3 -component  position  vector  of  the  center  of  mass  of  the  extendable  part  of  the  ith 
strut  defined  relative  to  the  {W}  origin;  i=l  ,6. 

(r  J 3-component  position  vector  of  the  ith  moving  platform-strut  joint  location  (Ps) 

defined  relative  to  the  platform’s  center  of  mass  (C),  specified  in  {W}  coordinates; 
i=l,6  [units:  (m,  m,  m)T] 

w Rp  3x3  orthonormal  rotational  matrix  defining  the  orientation  of  frame  {P}  relative 

to  frame  {W}  coordinates,  note  {W}  and  {F}  are  parallel  frames  for  all  time. 

' R 3x3  orthonormal  rotational  matrix  defining  the  orientation  of  frame  {P}  relative 

to  frame  {F}  coordinates. 
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w R,,  3x3  orthonormal  rotational  matrix  defining  the  orientation  of  frame  {F}  relative  to 
frame  {W}  coordinates. 

{s}(.  orthogonal  coordinate  system  {s2,s3,s23}.  affixed  to  strut  i with  origin  at  joint  bj. 

{s2,s3,sJ3}.  are  specified  in  {W}  coordinate  system. 

(s, ).  3-component  unit  vector  defining  the  first  joint  axis  of  the  Hooke  joint,  parallel  to 

x-axis  of  {W}  coordinate  system  specified  in  {W}  coordinates;  i=l,6;  total  of  six 
( s,  )i  unit  vectors,  one  for  each  strut. 

(s2  )i  3-component  unit  vector  orthogonal  to  (s3)j  and  (s, ),  vector  specified  in  {W} 
coordinates;  i=l,6;  total  of  six  (s2  )j  unit  vectors,  one  for  each  strut. 

( s 23 )i  3 -component  unit  vector  orthogonal  to  ( s 2 )j  and  ( s3  )j  vector,  direction  defined  by 
((s2).  x (s3).)  specified  in  {W}  coordinates;  i=l,6;  total  of  six  (s23)j  unit  vectors, 
one  for  each  strut. 

(s3).  3-component  unit  vector,  aligned  along  the  strut’s  longitudinal  axis  specified  in 
{W}  coordinates;  i=l,6;  total  of  six  (s3)j  unit  vectors,  one  for  each  strut. 

( y transpose  of  the  ( ). 

(r  ) (Ta  ) torsional  friction  in  strut-platform  joints  (Pi  and  bj),  moving  platform  and 

base  platform,  respectively;  direction  of  corresponding  torque  vector  is  aligned 
with  relative  angular  velocity  vector  in  the  joint. 

U 3x3  matrix  transformation  from  the  Euler  angular  rates  to  the  inertia  frame 
definition  of  the  platform’s  absolute  angular  velocity 
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(v5  );  absolute  velocity  of  the  center  of  mass  of  the  nonextendable  part  of  strut  in  inertial 
frame  {W}  coordinates. 

v absolute  velocity  of  the  center  of  the  mass  of  the  moving  platform  in  inertial  frame 
{W}  coordinates. 

(vE).  absolute  velocity  of  the  center  of  mass  of  the  extendable  part  of  strut  in  inertial 
frame  {W}  coordinates. 

{W}  coordinate  system  (Xw,  Yw,  Zw ):  a fixed,  Cartesian  inertial  reference  frame,  frame 
of  reference  for  all  absolute  vectors 

wk  the  parameter  set  to  be  identified  at  the  kth  iteration  step  in  its  estimation  (chapter  3) 

x position  of  the  cart  in  horizontal  plane  (Chapter  3) 

x a velocity  of  the  cart  in  horizontal  plane  (Chapter  3) 

yk  updated  output  (calculated  from  measured  states  (xk)  and  estimated  parameters 
( wk ))  for  kth  iteration  step 

z-x-z  Euler  platform  frame  {P}'s  absolute  orientation  is  defined  relative  to  frame  {F} 
by  ordered  z-x-z  Euler  anglular  rotations:  First  <j>  rotation  about  frame  {F}  axis  Z, 
then  9 rotation  about  new  X axis,  and  finally  rotation  \\i  about  new  Z axis  (units: 
radians) 

(aB)7  absolute  angular  acceleration  of  body  ‘B’  of  strut  defined  in  inertial  frame  {W} 
coordinates,  (aj  =(«,),  since  (a>J,  =(©,),. 

a c absolute  angular  acceleration  of  the  moving  platform  defined  in  inertial  frame  {W} 
coordinates 

(aJ  absolute  angular  acceleration  of  body  ‘E’  of  strut  defined  in  inertial  frame  {W} 
coordinates,  (aa).  = (a£),  since  (coj.  = (to  E). . 
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P 3-component  vector  consisting  of  the  z-x-z  Euler  angles,  defines  the  moving 
platform’s  absolute  orientation. 

P 3 -component  vector  consisting  of  the  time  rate  of  change  of  p , Euler  angular  rates 

for  the  moving  platform. 

P 3 -component  vector  consisting  of  the  time  rate  of  change  of  p . 

0X  absolute  angular  position  of  the  bottom  pendulum  (link  1)  (Chapter  3). 

6X  absolute  angular  velocity  of  the  bottom  pendulum  (link  1)  (Chapter  3). 

02  absolute  angular  position  of  the  top  pendulum  (link  2)  (Chapter  3). 

d2  absolute  angular  velocity  of  the  top  pendulum  (link  2)  (Chapter  3). 

(<*0,  absolute  angular  velocity  of  body  ‘B’  of  strut  defined  in  inertial  frame  {W} 
coordinates.  Note,  ((oE).  = (coB).  = (to  ).  for  link  i expressed  in  frame  {W} 

coordinates,  angular  velocities  of  bodies  Bj  and  Ej  which  are  constrained  to  rotate 
together  due  to  prismatic  joint  connecting  the  two  bodies. 

(!)  . = [u]p  the  absolute  platform  angular  velocity  defined  in  inertial  frame  { W} 
coordinates. 

(°o,  absolute  angular  velocity  of  body  ‘E’  of  strut  defined  in  inertial  frame  {W} 
coordinates.  Note,  (co;.  ).  = (co^ ).  = (to  ),  for  link  i expressed  in  frame  {W} 

coordinates,  angular  velocities  of  bodies  Bj  and  Ej  which  are  constrained  to  rotate 
together  due  to  prismatic  joint  connecting  the  two  bodies. 

(cos.)  absolute  strut  angular  velocity  defined  in  inertial  frame  {W}, 

)» = K 1 & )' + K 1 & )« + K. ),  (»» ), 
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where 


[cos  j angular  component  of  strut  i about  (s2);  axis 

(cos  JJ . angular  component  of  strut  i about  (s3).  axis 
(a)  ) angular  component  of  strut  i about  (s  ) axis 

\ *’.v23  //  V 23  // 

M,  absolute  strut  angular  velocity  defined  in  inertial  frame  {W}  with  (a>  I 

V ss,  // 
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CHAPTER  1 
INTRODUCTION 

1.1  History  of  Parallel  Kinematic  Machines 

Although  due  in  part  to  their  anthropomorphic  nature,  serial  manipulators 
comprised  of  open  chain  linkages  are  commonly  used  in  industry.  However,  since  the 
science  of  robotics  has  developed,  interest  in  parallel  kinematic  mechanisms  has  grown 
so  that  many  of  the  serial  robots’  actuation  transmissions  are  actually  comprised  of  closed 
chain  linkages  (the  general  class  in  which  parallel  kinematic  mechanisms  belong).  The 
Stewart  platform  [ 1 ] was  one  of  the  first  implementations  of  the  parallel  kinematic 
mechanism  and  has  been  widely  adapted  as  an  aircraft  simulator  (see  Figure  1.1). 
Although  it  has  been  claimed  by  researchers  that  parallel  kinematic 
mechanisms/machines  (PKMs)  have  higher  structural  stiffness,  and  proportionally 
distribute  the  payload  to  the  links,  these  attributes  are  highly  dependent  on  their 
configuration  design  (e.g.,  Stewart-Gough,  Delta-like,  Triaglide,  Linapod-like). 

Compared  to  serial  manipulators,  PKMs  have  the  advantage  of  relatively  simple  inverse 
kinematic  calculations  yielding  unique  solutions,  but  have  the  disadvantage  of  the 
forward  kinematics  resulting  in  high  order  polynomials  with  multiple  solutions. 

Using  application  driven  design  methods,  PKMs  can  be  realized  that  exhibit 
sufficient  structural  stiffness  and  high  speed  motions  for  a wide  variety  of  applications 
(e.g.,  motion  simulations,  assembly  and  packaging  end-effectors,  surgical,  MEMS, 
entertainment).  Limited  usage  of  PKMs  is  largely  due  to  small  workspace,  small  regions 
of  dexterity,  variable  end-effectors  stiffness,  non-unique  solutions  of  forward  kinematics. 


1 


2 


workspace  singularities,  self-motions  and  uncomprehended  dynamic  performances. 
During  the  last  decade,  there  have  been  many  debates  on  the  ‘true’  attributes  of  PKMs 
and  their  applications,  especially  since  their  1 994  debut  as  a new  machine  tool  paradigm. 
In  the  literature  [2],  it  was  described  as  follows. 

The  PKM  suffers  from  a lack  of  comprehensive  calibration  techniques,  from  weak 
or  not  well-defined  dynamics  behavior,  and  from  the  need  for  complex  control  algorithms 
for  the  compensation  of  stiffness  variations  in  the  workspace. 

To  address  the  above  issues  and  subsequently  enable  greater  understanding  and 
utilization,  this  dissertation’s  research  focuses  on  the  dynamics.  The  next  few  sections 
include  the  brief  summary  of  the  research  accomplishments  and  applications,  which  are 
now  being  implemented,  followed  by  the  proposed  methodology,  and  research  plan. 

1.2  Types  of  PKM  and  Applications 

Parallel  mechanisms  can  be  classified  as  3 basic  types: 

• Variable  link  length  (Stewart-Gough) 

• Variable  two  link  orientation  (Delta-like) 

• Fixed  link  length  with  variable  base  joint  positions  (Triaglide  or  Linapod- 
like) 

These  types  of  PKMs  and  their  variants  are  commonly  found  in  industry  and  research 
laboratories.  The  following  highlights  a few  of  the  variants  in  design. 

For  a spherical  parallel  manipulator  (SPM)  designed  to  orient  a rigid  body  in  space, 
researchers  used  kinematic  conditioning  as  a performance  index  to  optimize  the 
manipulator  configuration  [3,  4],  Figure  1.2  shows  a sketch  of  the  SPM.  This  manipulator 
could  be  used  for  the  orientation  of  parabolic  antennas,  panels,  and  machine  tool  beds. 
Such  a mechanism  can  also  be  applied  as  an  orientation  of  the  wrist  for  an  industrial 
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robot.  The  research  by  Gosselin  and  Angles  [3]  includes  the  restriction  of  symmetry  of 
the  machine  kinematics,  workspace  maximization,  and  isotropy  of  the  machine.  A 
discussion  of  the  singularities  was  also  included. 

Many  researchers  proposed  novel  kinematic  structures  for  a parallel  manipulator 
with  6 DOF  (Figure  1.3).  Byun  and  Cho  [5]  showed  a feasible  real-time  computation  for 
all  forward  solutions  of  a 6 DOF  parallel  manipulator  consisting  of  a base,  platform,  and 
3-PPSP  (P-Prismatic,  S-Spherical)  serial  mechanisms.  Their  machine  was  to  overcome 
the  common  drawbacks  of  the  conventional  six-degree  of  freedom  parallel  manipulators, 
which  are  complex  forward  computation,  small  workspace  volume,  and  limited  load 
capacity.  However,  the  resulting  PKMs  are  still  constrained  by  a relatively  small 
functional  workspace.  These  PKM  structures  could  be  used  as  a motion  base,  which 
generally  needs  only  a small  workspace  capacity.  The  authors  [5]  did  not  address 
workspace  restrictions,  singularity  points,  and  the  dynamics  of  the  structure. 

Other  researchers  have  conceived  an  idea  for  a parallel  manipulator  (Turin  parallel 
manipulator  named  TuPaMan,  see  Figure  1 .4)  [6],  This  manipulator  is  composed  of  three 
articulated  double  parallelograms,  three  prismatic  joints,  three  spherical  joints,  and  a 
moving  upper  platform.  The  TuPaMan  was  strongly  influenced  by  the  Stewart  platform 
and  can  be  used  as  a manipulator  wrist  or  arm.  The  restrictions  for  TuPaMan  are  a 
function  of  the  size  of  the  device,  using  the  overall  dimensions  and  the  area  of  the  middle 
cross-section  of  the  workspace  as  synthetic  performance  indices.  Oh  [7]  used  linear 
motors  to  extend  the  workspace.  He  presents  a method  to  overcome  three  constraints: 
strut  collision,  base  puck  collision,  and  singularity  considerations. 


4 


1.3  PKM  Kinematics 

The  PKM  kinematics  has  been  extensively  investigated  for  the  past  few  decades. 
Much  of  kinematic  research  focus  has  been  concentrated  on  solving  the  forward 
kinematics  problem.  The  traditional  Stewart  Platform  Mechanism  has  40  mathematically 
achievable  configurations  in  the  solution  set.  Duffy  et  al.  [8]  reduced  the  forward 
kinematic  solutions  to  eight  for  a special  6-6  PKM. 

1.4  PKM  Dynamics  and  Controls 

To  date,  not  much  research  on  the  dynamics  of  the  PKM  has  been  done.  The 
dynamic  equations  of  motion  have  been  derived  for  the  generalized  Stewart  platforms 
using  Lagrangian,  Newton-Euler,  and  Kane’s  method.  Of  the  work  available  in  the 
literature,  a sampling  of  researchers  who  investigated  modeling,  dynamics  and/or 
parameter  identification  in  the  robotics  community  are:  [9-23].  Sestieri  [9]  derived  a 
series  of  simple  linear  problems  from  the  non-linear  problem,  and  solved  them  with  a 
frequency-response  function  (FRF).  Tinone  and  Aoshima  [10]  used  a repetitive  controller 
in  order  to  reduce  bias  from  measurement  error.  Lin  [16]  used  MLC  (minimal  linear 
combinations)  for  avoiding  the  difficulties  to  estimate  the  parameters.  The  MLC  was 
used  for  estimating  the  parameters  in  the  serial  manipulators,  but  it  can  also  be  used  in 
the  parallel  manipulator.  Lin  described  MLC  excluding  the  friction  parameters.  Further 
insight  into  PKM  and  serial  robot  dynamics  and  identification  has  been  provided  in  the 
following:  vibration  analysis  [24],  trajectory  analysis  [25],  force  analysis  [26]  accuracy 
analysis  [27],  forward  analysis  [28],  stiffness  analysis  [29],  and  parameter  estimation  [BO- 
SS]. Available  in  the  literature,  only  one  reference  [33]  attempted  parameter 
identification  of  PKM’s  rigid  body  dynamics. 
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1.5  Prior  Work  within  SAMM  Lab 

The  following  is  a summary  of  the  research  on  PKMs  conducted  within  the  Space, 
Automation  and  Manufacturing  Mechanisms  (SAMM)  Laboratory.  Hardage  [34] 
developed  a methodology  for  identifying  the  structural  dynamic  parameters  of  PKMs 
validated  using  experimental  data  from  the  Hexel  PKM  located  at  Sandia  National 
Laboratory.  He  used  component  mode  synthesis  and  Lagrangian  techniques  to  derive  the 
equations  of  vibratory  motion.  Modal  testing  was  performed  to  measure  the  stiffness  and 
frequency  of  the  Sandia  National  Laboratories  Hexel  PKM.  His  results  demonstrated  that 
the  Hexel  machine  could  not  be  expected  to  achieve  high  structural  performance  nor  be 
widely  adopted  as  a machine  tool. 

Johnston’s  research  [35]  was  to  develop  dynamic  performance  indices  for  a parallel 
manipulator.  In  his  work.  Kane’s  method  was  used  to  establish  the  equations  of  dynamic 
motion.  To  eliminate  mixed  units  in  this  metric,  he  modified  the  energy  expression  to  be 
in  terms  of  linear  velocities  of  the  platform  only.  The  performance  indices  clearly 
demonstrated  that  to  fully  characterize  the  dynamic  performance  of  the  machine, 
kinematics  based  indices  are  not  sufficient.  This  was  further  substantiated  by  Walker’s 
research  on  studying  the  effects  of  kinematic  errors  on  dynamic  performance. 

Register  [36]  characterized  the  effects  of  friction  using  a parameter  estimation 
algorithm  (extended  Kalman  filter).  In  his  work,  friction  forces  were  shown  to  be 
dominant  over  inertial  and  gravitational  forces  in  linear  and  angular  motions  of  the  struts. 
Estimating  the  friction  parameter  for  the  single  joint  strut  was  achieved  within  10  percent 
using  the  estimation  algorithm,  however,  the  identification  of  friction  parameters  in 
multiple  joint  types  was  not  addressed  in  his  thesis.  Scott  Shamblin  [37]  conducted 
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further  analysis  and  experimental  identification  of  the  individual  struts  and  the  assembled 
PKM  and  again  was  faced  with  the  challenges  of  high  friction. 

1.6  System  Identification  Methods 

System  identification  methods  are  well  known  and  utilized  by  researchers  working 
in  the  areas  of  dynamics,  controls,  and  estimation  theory.  System  identification  can  be 
categorized  into  state,  parameter  and  dual  estimation.  The  parameter  estimation  method  is 
used  for  approximating  the  unknown  parameters  of  the  system  from  measured  data  which 
contains  the  known  system  input  and  output  variables.  There  are  many  methods  existing 
in  parameter  estimation,  however  due  to  its  mathematical  simplicity,  a commonly  used 
method  is  least  squares.  The  least  squares  method  is  widely  used  to  find  system 
parameter  and  state  estimation.  Least  squares  methods,  however,  exhibits  convergence 
problems  in  the  presence  of  noise  and  lacks  the  ability  to  effectively  handle  nonlinearities 
in  system  identification.  This  has  lead  to  the  development  of  various  mutations  of  least 
squares  [23],  Of  these  mutations,  a recent  method  gaining  popularity  is  the  extended 
Kalman  filter.  The  most  efficient  way  to  identify  nonlinear  system  parameters,  however, 
is  the  unscented  Kalman  filter  method  (UKF)  [38],  For  this  dissertation’s  research,  the 
unscented  Kalman  filter  method  is  selected  as  the  method  of  choice  for  dealing  with  the 
dynamic  complexity  of  PKMs.  The  major  advantage  for  the  UKF  is  that  it  does  not 
require  the  calculation  of  Jacobian  and  Hessian  matrices.  This  approach  will  be  shown  in 
Chapter  4 to  yield  faster  and  more  convenient  solutions  for  the  nonlinear  equations  of  the 
PKM.  An  example  for  comparing  the  least  squares  method  and  the  unscented  Kalman 
filter  method  is  provided  in  Chapter  3,  demonstrating  the  effectiveness  of  the  UKF 


method. 
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1.7  Problem  Statement 

When  a system  has  unknown  inertial  parameters  that  are  difficult  to  identify  by 
simple  physical  measurements  of  the  components,  a parameter  estimation  method  can 
resolve  the  variables  using  commanded  inputs  and  measured  output  response  values.  The 
main  goal  of  this  dissertation  is  to  develop  a methodology  for  identifying  the  system 
dynamic  model  of  PKMs  through  estimating  the  inertial  parameters.  The  important 
characteristic  of  the  inertial  parameters  of  the  machine,  in  general,  is  that  those  values  are 
invariant  in  time.  In  other  words,  when  changing  the  state  values  of  system  over  time, 
the  inertial  parameters  stay  the  same  regardless  of  what  the  input  and  output  values. 

The  outline  of  the  procedure  established  to  obtain  the  inertial  parameters  presented 
in  this  dissertation  is  as  follows: 

First,  formulate  a model  base  by  deriving  the  equations  of  motion.  The  theoretical 
background  relating  to  the  PKM  dynamics  and  estimation  method  is  shown  (Chapter  2). 
Second,  analyze  and  compare  candidate  estimation  techniques.  Least  squares  (LS) 
method,  extended  Kalman  filter  (EKF),  and  unscented  Kalman  filter  (UKF)  were  selected 
for  comparisons  as  representing  the  spectrum  of  identification  methods  (Chapter  3).  For 
this  comparison,  the  inverted  double  pendulum  example  was  selected.  Simulation  results 
were  also  generated  for  exploring  both  force-based  modeling  and  energy-based  modeling 
as  a model  basis  for  implementing  the  preferred  estimation  method  (UKF,  Chapter  4). 
Third,  experimental  validation  of  the  capability  and  the  implementation  of  the  UKF 
estimation  method  is  demonstrated  in  Chapter  5.  The  measured  PKM  input  and  output 
responses  from  the  experiment  were  employed  in  the  estimation  program.  Finally,  the 


validation  and  reconciliation  of  the  resulting  parameter  identification  methodology  were 
completed. 
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Figure  1.3  Novel  6-DOF  Parallel  Figure  1.4  Turin  Parallel  Manipulator 

Machine  [5]  (TuPaMan)  [6] 
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Figure  1.5  Generalized  Stewart  Platform  Schematic  [36] 


CHAPTER  2 

THEORETICAL  BACKGROUND 


Each  of  the  parameter  identification  methods  considered  in  this  dissertation  require 
as  a model-base  the  theoretical  model  of  the  expected  system  dynamics.  In  the  derivation 
of  the  theoretical  model,  two  different  approaches  were  taken  in  formulating  a model- 
base.  One  approach  uses  the  system’s  equations  of  motion,  referred  to  as  the  force-based 
model.  The  other  approach  uses  the  expressions  from  the  system’s  energy  equations, 
referred  to  as  the  energy-based  model.  The  same  force-based  model  also  serves  as  a 
theoretical  representation  of  the  physical  system  in  the  simulated  identification 
implementation  studies  and  other  supporting  studies.  In  this  chapter,  the  theoretical 
modeling  of  the  general  class  of  Stewart-Gough  parallel  kinematic  mechanisms  (PKM)  is 
presented  with  supplemented  details  given  in  Appendices  A and  B.  The  chapter 
continues  with  presentation  of  the  theoretical  details  of  three  parameter  identification 
methods.  Simulated  and  experimental  implementations  of  these  methods  are  presented  in 
Chapters  3,  4 and  5.  The  chapter  closes  with  a discussion  of  implementation  issues  and 
preliminary  characterization  of  the  expected  dynamic  parameter  dominance,  the  ease  in 
identifying  various  terms,  and  effects  of  limited  and/or  indirect  sensing.  While  the 
dissertation  results  are  limited  to  the  Stewart-Gough  class  of  PKMs,  it  should  be  noted 
that  the  methodology  developed  in  this  dissertation  and  the  issues  addressed  are  directly 
applicable  to  PKMs,  in  general. 


12 


2.1  System  Description  and  Model 

In  this  dissertation,  the  Stewart-Gough  PKM  type  is  assumed.  This  type  of  PKM 
has  six  extendable  struts  (linear  actuators)  that  provide  6-degrees  of  freedom  (d.o.f.)  to  a 
moving  platform  relative  to  a base  platform,  refer  to  schematic  in  Figure  2. 1 . Figure  2.2 
shows  a physical  diagram  of  a representative  strut,  and  figure  2.3  is  its  corresponding 
lumped  mass  model  where  details  are  provided  in  the  section  on  equations  of  motion. 

The  following  describes  three  PKM’s  that  are  a part  of  the  Stewart-Gough  family  of 
PKM’s.  Each  consists  of  different  mass  distribution  and  special  joint  arrangements  where 
the  strut  design  is  generally  chosen  by  the  designer  to  minimize  computation  and  to 
facilitate  actuation.  In  section  2.4,  these  three  PKM  configuration  variations  were  used  to 
characterize  the  expected  parameter  dynamics  and  to  delineate  potential  implementation 
issues. 

2.1.1  System  Description 

2.1.1. 1 UF-PKM 

Figure  2.4  shows  the  Special  6-6  PKM  [8]  which  is  the  design  type  of  the 
University  of  Florida  - Space,  Automation  and  Manufacturing  Mechanisms  (UF-SAMM) 
Laboratory’s  testbed’s  PKM,  a subset  of  Stewart-Gough  PKMs.  The  special  6-6  PKM 
consists  of  a moving  platform  with  its  6 degrees  of  motion  actuated  via  six  extendable 
struts.  The  special  6-6  PKM  has  a unique  telescoping  strut  arrangement  such  that  both 
the  platform  and  base  consist  of  three  sets  of  three  joints  lying  on  lines  that  form  the 
sides  of  a triangle.  It  is  this  joint  configuration  that  yields  reduced  computation 
requirement  for  the  forward  kinematic  solutions  compared  to  that  required  of  the  general 
Stewart-Gough  PKM  as  well  as  allows  for  simpler  physical  joint  designs.  As  shown  in 
Figure  2.2,  each  strut  consists  of  a motor,  encoder,  brake,  spherical  joint,  ball  screw, 
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Hooke  joint  (also  known  as  universal  joint  or  cardan  joint),  hollow  ball  spline/nut,  and 
outer  tube.  Each  of  the  struts  (linear  actuators)  is  connected  to  the  stationary  base 
platform  through  a large  hollow  spherical  joint.  This  connection  allows  the  strut  to  rotate 
through  space  as  its  length  is  varied.  A servomotor,  whose  armature  is  connected  to  a 
ball  screw  that  passes  through  the  center  of  the  hollow  spherical  joint,  drives  the  strut.  To 
prevent  rotation  about  the  servomotor  axis,  the  combination  of  a hollow  ball  spline  and  a 
Hooke  joint  is  used  as  the  connection  between  the  moving  platform  and  the  strut. 

2. 1.1.2  SNL-PKM 

Figure  2.4  (b)  presents  the  Sandia  National  Laboratory’s  Hexel  parallel  machine 
tool.  Still  a Stewart-Gough  type  PKM,  its  configuration  differs  in  that  the  motor  and  ball 
nut  of  each  strut  are  concentric  with  the  base  spherical  joint  and  the  ball  screw  passes 
through  the  base  joint  connecting  directly  to  the  moving  platform,  thus  eliminating  the 
need  for  a hollow  ball  spline  and  outer  tube.  The  Hooke  joints  are  replaced  with 
bifurcated  ball  joints  that  allow  for  concentric  spherical  joint  motion  between  two  struts 
and  the  moving  platform.  The  joints  are  arranged  between  the  struts  and  the  base  and 
moving  platforms  as  illustrated  in  Figure  2.4  (b).  Between  the  moving  platform  and 
struts,  the  joint  arrangements  form  an  equilateral  triangle. 

2.1. 1.3  FORD-PKM 

Figure  2.4  (c)  is  of  the  Ford  PKM  used  for  full  6-degrees  of  freedom  vehicle 
simulation  with  the  driver  in-the-loop  control.  The  Ford  PKM  has  an  approximate 
octahedral  layout  similar  to  the  Hexel  base  platform.  The  platform  joint  configuration(s) 
are  arranged  in  close  pairs  for  both  the  base  and  moving  platforms.  This  PKM  is 
powered  by  hydraulics,  i.e.,  each  strut  is  essentially  a hydraulic  cylinder  with  its  ends 
connected  to  the  base  and  moving  platforms.  Both  the  Ford  and  UF  PKMs  have 
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relatively  short  struts  in  relation  to  their  base  joint  spread,  but  are  totally  different  in 
terms  of  the  mass  distribution.  I.e.,  the  moving  platform  of  the  UF  PKM  has  relative 
small  inertia  compared  to  the  inertia  of  struts,  while  the  Ford  PKM  has  the  majority  of  its 
inertia  in  the  moving  platform. 

2.1.2  Lumped  Mass  Parameter  Model 

The  physical  system  has  three  major  components  to  be  categorized:  the  stationary 
base,  the  strut  and  the  moving  platform.  The  stationary  base  serves  as  the  inertial 
reference  frame  for  the  system.  The  moving  platform  has  six  degrees  of  freedom  and 
consists  of  a single  mass.  It  has  a mass  mP  and  inertia  tensor  IP.  The  last  major 
component  is  the  strut.  Since  the  actuators  only  produce  translational  motion  along  the 
longitudinal  axis  of  the  struts,  the  actuating  component  is  modeled  as  a linear  force 
applied  across  a prismatic  joint.  For  the  UF-PKM,  the  struts  are  connected  between  the 
stationary  base  platform  and  a moving  platform  with  a spherical  joint  and  a Hooke  joint, 
respectively.  However,  without  loss  of  generality  and  for  modeling  purposes,  each  strut 
can  be  theoretical  modeled  as  forming  a kinematic  chain  (H-P-S)  which  consists  of  a 
Hooke  joint,  a prismatic  joint  (representing  the  linear  actuation),  and  a spherical  joint. 
Both  kinematic  chains  have  the  same  kinematic  degrees  of  freedom.  Figures  2.2  and  2.3 
show  the  physical  diagram  of  a single  strut  for  the  UF  PKM  and  the  corresponding 
lumped  mass  model.  For  Body  Bj,  the  (mb)i  and  (lb),  are  respectively  the  composite  mass 
and  inertia  tensor  for  all  the  components  in  the  physical  model  that  do  not  extend  (or 
translate  along  strut  axis  (s3)0  as  the  strut  (linear  actuator)  changes  in  length.  For  the 

UF-PKM,  this  consists  of  an  outer  tube,  ball  screw,  ball  of  the  spherical  joint,  motor, 
encoder  and  brake.  The  strut  length  is  defined  by  the  distance  between  center  of  the 
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moving  platform-strut  joint  (P;)  and  base-platform-strut  joint  (bj).  For  Body  Ei?  the  (me)i 
and  (Ie)i  are  respectively  the  composite  mass  and  inertia  for  the  remaining  components 
that  exhibit  a translational  motion  along  the  ( s3  )i  axis.  For  the  UF-PKM,  this  consists  of 

a ball  nut,  hollow  ball  spline,  and  Hooke  joint.  The  actuator’s  force  (Fa)j  can  be 

modeled  as  a force  exerted  between  these  two  lumped  mass  entities.  The  joint  friction 
within  the  kinematic  chain  can  be  modeled  as  both  translational  and  torsional  friction, 

(F(i  )i,  (7^  )i,  and  ( 7^  )b  respectively.  The  quantities  ( Ln );,  (LE)h  and  Lj  denote 
respectively  the  position  along  strut  axis  (s3)j  relative  to  point  (bj)  of  the  lumped  masses, 
(me)j  and  (mb)j,  center  of  mass  and  Li  is  the  instantaneous  length  of  strut  i.  In  reiterating 
the  above  and  illustrated  in  Figure  2.3,  the  strut  masses  are  lumped  into  two  parts.  Body 
E j (me,  Ie)  j and  Body  B , (mb,  lb)  j,  which  both  experience  rotational  motion  relative  to  the 
base  joint,  i.e.,  rotate  with  (s3)j  axis.  In  addition,  Body  Ej  experiences  translational 

motion  along  (s3)j  where  Body  Bj  does  not.  Since  Body  Ej  translates  freely  in  the 
direction  of  the  strut  axis,  the  location  of  its  center  of  mass  is  located  at  a variable 
distance  measured  from  the  base  platform-strut  joint  location  (bj)  along  the  strut  axis 

(s3  )j.  As  stated  earlier,  this  distance  is  denoted  by  ( LK  )\.  The  velocity  and  acceleration 

of  Body  Ej  can  be  expressed  in  terms  of  the  linear  velocity  and  acceleration  of  the 
extendable  part  of  the  strut  in  the  direction  of  the  strut  axis  and  the  strut  angular  velocity 
and  acceleration  about  axes  defined  in  the  base  joint.  Body  Bj  has  a center  of  mass  that  is 

located  a fixed  distance  ( LB)\  from  the  base  joint  (bj)  and  is  measured  in  the  direction  of 
the  strut  axis.  Since  this  distance  is  fixed,  the  velocity  and  acceleration  of  Body  Bj  can  be 
expressed  in  terms  of  the  strut’s  angular  velocity  and  acceleration  about  the  axes  defined 
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in  the  base  joint  (bj).  To  summarize,  the  total  number  of  the  lumped  masses  in  the 
Stewart-Gough  PKM  model  are  thirteen:  two  for  each  strut  and  one  for  the  moving 
platform. 

2.1.3  Coordinate  System  Assignment 

In  the  analysis  of  kinematics  and  dynamics,  it  is  typical  that  a number  of  coordinate 
systems  are  defined  to  describe  the  motion  of  bodies.  To  formulate  the  dynamic 
equations  in  this  dissertation,  three  types  of  coordinate  systems  are  defined:  global  (the 
inertial  reference  frame),  strut  (S-coordinate  system),  and  moving  platform  coordinate 
systems  (translational  and  rotational  coordinate  system). 

2. 1.3.1  Global  and  Moving  Platform  Coordinate  Systems 

The  global  (inertial)  coordinate  system  is  denoted  as  {W}  coordinate  system.  It  is  a 
fixed,  Cartesian  inertial  reference  frame.  This  coordinate  system  can  be  located 
arbitrarily,  however,  for  the  PKMs  considered,  the  {W}  coordinate  system  is  fixed  at  the 
geometric  center  location  of  the  stationary  base  platform  as  shown  in  Figure  2.5.  The 
moving  platform  coordinate  system  consists  of  two  frames:  a translational  and  rotational 
coordinate  system.  The  translation  frame  is  called  the  {F}  coordinate  system  and  the 
rotation  frame  is  called  the  {P}  coordinate  system.  As  in  Figure  2.5,  the  origin  of  the  {F} 
coordinate  system  moves  with  the  center  of  mass  location  of  the  moving  platform  but  the 
orientation  of  the  {F}  frame  remains  parallel  with  the  { W}  frame  for  all  time.  It  shares 
the  same  origin  with  the  {P}  coordinate  system.  The  {P}  frame  is  affixed  to  the  moving 
platform,  where  its  absolute  orientation  and  position  is  used  to  describe  the  motion  of  the 
moving  platform  relative  to  the  {W}  frame.  The  {F}  frame  is  used  as  an  intermediate 
step  for  obtaining  the  {P}  frame’s  absolute  position  and  orientation.  The  {P}  frame 
orientation  is  arbitrary,  but  is  typically  selected  to  have  its  axes  aligned  with  the  principal 
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axes  of  inertia  (as  done  in  this  dissertation’s  examples).  To  avoid  singularities  in  the 
computation  the  axes  are  also  selected  to  have  specific  orientations  when  the  platform  is 
in  its  ‘home’  position.  This  is  explained  below  in  section  2. 1 . 1 .2  c. 

The  naming  convention  used  for  denoting  the  center  location  of  the  strut  joint  on 
the  base  and  moving  platform  as  well  as  the  center  of  mass  location  of  the  moving 
platform  are  shown  in  Figure  2.6.  All  joints  on  the  stationary  base  are  labeled  point  bj 
where  i=l,6  for  six  struts.  Their  locations  are  defined  in  the  {W}  coordinates  (e.g.,  wb  )• 

On  the  moving  platform,  the  locations  of  the  joints  are  denoted  as  Pi.  The  displacement 
from  the  origin  of  the  {F}  or  {P}  coordinate  system  to  point  Pj  is  called  vector  (wr  F, 

i.e.,  location  of  Pi  relative  to  the  platform’s  center  of  mass,  C,  expressed  in  {W} 
coordinates. 

2. 1.3.2  Strut  Coordinate  System 

The  strut  coordinate  system  is  denoted  as  the  { S } i coordinate  system  (refer  to 
Figure  2.7  which  shows  the  theoretical  Hooke  joint  model  for  the  base  joint  bj).  The  (s,  )j 
unit  vector  is  used  to  denote  the  first  axis  in  the  Hooke  joint.  The  angular  rate  about  this 
axis  is  defined  as  (CO  )j.  Used  as  a reference  to  define  other  vectors  and  without  loss  of 

sl 


generality,  the  direction  of  the  (s,  )i  vector  is  arbitrarily  defined  to  be  parallel  to  the  Xw- 
direction  of  {W}  coordinate  system.  I.e.,  the  X-axis  of  the  X-Y-Z  coordinate  frame  in 
Figure  2.7  theoretically  remains  parallel  to  the  {W}  coordinate  system  for  all  time.  The 
(s3 ) unit  vector  is  defined  as  the  unit  vector  along  the  line  of  action  between  joint  Pi  and 
bj,  which  is  calculated  as 


(U 


CU-fO, 

!('<•,), -('0,1 


(2-1) 
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Note,  the  subscripts  i denote  ith  strut  where  i=l,6  for  the  PKM’s  under  consideration. 

The  (s  2 ).  unit  vector  is  formed  as  the  cross  product  of(s, ) and  (st ) divide  with  the 

magnitude  of  the  cross  product  of  (§, ).  and  (s3 ) . The  (s23 ) unit  vector  is  defined  as  the 
cross  product  of  (s2 ).  and  (s3 ) . Thus,  the  three  orthogonal  unit  vectors  form  the  strut 


coordinate  system  expressed  by  equations  (2.1)  through  (2.3)  in  {W}  frame  coordinates. 


(s.),=[l  0 Of  and  (g  ) =fei-XS4 
(S»),  =(S.XS,)( 


S,  X slU 

I 1 J 11/ 


(2.2) 

(2.3) 


where  the  strut  coordinate  system  {s}.  = |s2,s3,s23 } . - This  frame  definition  is  used  for 

describing  the  strut  motion  and  to  facilitate  the  derivation  of  the  strut  dynamic  terms  in 
the  equations  of  motion. 


2.2  Equations  of  Motion 

To  conduct  model-based  dynamic  parameter  identification,  one  must  first  formulate 
the  equations  of  motion  used  to  represent  the  expected  dynamic  model  of  the  system. 

The  equations  used  in  this  dissertation  are  based  on  prior  theoretical  modeling  by 
references  [35,  36,  and  39-41],  Assuming  a lumped-mass  rigid  body  model,  Kane’s 
method  is  selected  for  deriving  the  dynamics  equations  of  motion.  The  following 
summarizes  the  assumptions  made. 

• Rigid  body  dynamics  of  a system  of  lumped  masses  and  inertias. 

• Strut  inertias  are  symmetric  about  strut  longitudinal  axis  defined  along  line  of 
action  between  the  strut  base  and  platform  joint  centers. 

• No  kinematic  joint  errors. 

• PKM  does  not  have  kinematic  singularities  within  its  workspace  (volume),  e.g., 
det[J]  and  det[U]  are  not  equal  to  zero  for  all  configurations  (these  terms  are 
defined  later  in  this  chapter). 
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• Center  of  mass  locations  of  the  extendable  portion  of  the  strut  (Body  E)  j and 
nonextendable  portion  of  the  strut  (Body  B)i  remain  along  (s3)j.  Both  (me)j  and 
(mb)i  rotate  with  (s3)j  axis  (Refer  to  Figure  2.2  and  2.3). 

• Unmodeled  dynamics  are  attributed  to  joint  imperfections,  mass  variations,  actuator 
and  controller  dynamics  and  the  inertial  effect  of  rotational  motion  of  struts  about 
the  struts  longitudinal  axis.  These  are  generally  assumed  negligible  but  are  retained 
in  the  theoretical  model  to  denote  sources  of  error  in  parameter  identification. 

• The  struts  are  designed  to  restrict  free  spin  rotation  about  its  longitudinal  axis  (s3)j. 

• Modified  joint  friction  model  ( F(i  )j  that  includes  both  viscous  and  Coulomb 

friction  is  assumed  [36]  (model  details  are  provided  in  Chapter  4).  In  the  analysis 
presented,  the  torsional  friction  within  the  joints  between  the  base  and  moving 
platforms  and  the  struts  are  assumed  negligible  (i.e.,  ^ j = (f.  ) = 0>  i=l ,6). 

2.2.1  Generalized  Coordinates 

To  describe  the  motion  of  the  6 degrees  of  freedom  PKM  system,  a set  of  6 
generalized  coordinates  must  be  chosen.  In  this  system,  the  generalized  coordinates  are 
chosen  to  be  the  center  of  mass  location  of  the  moving  platform  (xc,  yc,  zc)  with  reference 
to  the  { W}  coordinate  system.  For  the  moving  platform’s  absolute  orientation,  the  z-x-z 
Euler  angles  are  chosen  to  avoid  computational  difficulties  and  to  allow  direct  integration 
of  the  angular  velocity.  The  z-x-z  Euler  angles  represents,  sequentially,  a rotation  (j) 
about  the  z-axis  of  the  {W}  coordinate  system,  a rotation  0 about  the  x-axis  of  the  rotated 
coordinate  system,  and  a rotation  vp  about  the  z-axis  of  the  rotated  coordinate  system.  For 
this  system,  to  avoid  the  mathematical  singularity,  when  the  platform  orientation  is  in  its 
‘home’  position  (parallel  to  the  base  frame  {W}  and  at  a specified  (xc,  yc,  zc)),  the  z-x-z 
Euler  angles  are  set  to  ((j),  0,  vj/)  = (0°,  -90°'  0°).  Thus,  at  ‘home’  position,  the  {P}  frames 


are  expressed  as  {Yp=  -ZF,  Xp=  XF,  Zp=  YF}. 


The  symbol  q is  used  to  express  the  6x1  vector  of  the  generalized  coordinates.  For 
Kane’s  equations,  generalized  speeds  are  also  required.  The  generalized  speeds  ( q ) 
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correspond  to  the  time  derivative  of  the  generalized  coordinates.  The  generalized 
coordinates  and  the  generalized  speeds  can  be  expressed  as  follows,  respectively. 


q = k 

(2.4) 

q = [xc,yc,zc,0,$,yr] 

(2.5) 

where  q is  the  absolute  position  and  orientation  coordinates  of  the  platform  at  its  center 

of  mass  denoted  by  point  C in  Figure  2.6.  In  equation  (2.5),  the  ( ) denotes  the  time 
derivative  of  ( ).  The  generalized  coordinates  can  be  further  expressed  as  q = [w  r T , pT  j 
where  the  platform’s  absolute  position  (wr  ) referenced  in  {W}  is 

"r C=[xc  yc  zcJ  (2.6) 

The  orientation  coordinates  of  the  moving  platform  can  also  be  expressed  by  a vector  of 
the  Euler  angles,  p , where 

P = [^  & yX  (2.7) 

Similarly,  generalized  speed  can  be  expressed  as 

q = [v.'-P1  J (2.8) 

where  v is  the  absolute  linear  velocity  of  the  platform’s  center  of  mass. 

2.2.2  Kinematics 

In  order  to  obtain  the  dynamic  response  of  the  PKM  system,  the  position,  velocity, 
and  acceleration  of  all  bodies  must  be  expressed  in  terms  of  the  generalized  coordinates. 
However,  to  obtain  the  position,  velocity  and  acceleration  of  Bodies  Bj  and  E,  in  terms  of 
generalized  coordinates,  the  inverse  kinematics  must  be  solved.  The  inverse  kinematics 
is  that  given  position  and  orientation  data  of  the  moving  platform,  one  can  determine  the 
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length  ( Z( ) and  orientation  (s3 )(.  of  the  struts.  The  subscript  i in  these  and  the  following 

expressions  and  equations  denote  that  there  is  a corresponding  expression  and  set  of 
numerical  values  for  each  strut  i for  i=l,6.  Without  loss  of  generality,  the  i subscript  will 
be  dropped  where  it  is  understood  that  for  all  joint,  there  are  i=l,6  set  of  expressions. 
Given  the  orientation  of  the  moving  platform,  a standard  3x3  rotation  matrix  is  defined  to 
map  the  orientation  from  the  {P}  coordinate  system  to  {F}  or  {W}  coordinate  systems 
per  the  following. 


w 


C\]JC(f)  - C0S^V|/ 
CVj/5^  + C0C^\|/ 
j-V|/S0 


- si|/c^  - cB^^cii/ 

- si)/^  + c0c$!cy|/ 

Cl)/.?© 


S0S^ 
- sQc(f> 

C0 


(2.9) 


where  s and  c denote  the  sine  and  cosine  functions,  respectively,  and  " = [I]  since 

{F}  and  {W}  are  parallel  frames. 

Since  the  position  of  the  strut  joint  with  the  moving  platform  joint  (P,)  is  denoted 
by  (r pc ).  with  reference  to  frame  {P},  pre-multiplying  {rpc ) by  w Rp  yields  the  value  of  Pj 

position  (”  rp).  in  the  {F}  coordinate  system.  Then  to  complete  the  mapping  to  {W} 


frame,  the  translation  from  the  origin  of  {F}  frame  to  {W}  frame  (wr  ) must  be  done. 
The  total  mapping  of  Pj  is  shown  in  equation  (2.10). 


r ) = 


r +“R 


'(rpc), 


(2.10) 


Once  the  location  of  the  position  Pj  is  defined,  then  the  strut  length  and  orientation 
can  be  determined.  Subtracting  the  location  of  point  bi  from  the  location  of  Pi  gives  the 


length  of  the  strut. 

(Zs3),  =CrP-wrh)i 


(2.11) 
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Since  the  unit  vector  of  the  strut  ( s3  )j  expressed  in  { W}  coordinates  is  aligned  to 
the  individual  strut  i,  it  can  be  found  with  the  following  equation. 


r r, -*!*), 


(S3);  = 


w r _wr  I 

lP  '/> 


(2.12) 


The  velocity  of  point  Pj  ((v  )()  can  be  calculated  by  differentiating  (H  r;.)  in 
equation  (2.10)  with  respect  to  time. 


(v  ) = v + (to  x"R  r ) = v +(to  x"r  ). 

^ p ' 1 c ' c p pc  / / c \ c pc  / I 


(2.13) 


where  (Oc  is  the  moving  platform’s  absolute  angular  velocity  defined  in  the  inertial 
frame  {W}  coordinates. 

Since  the  velocity  of  point  P must  be  expressed  in  terms  of  the  generalized 
coordinates,  an  expression  for  (Or  in  terms  of  the  generalized  coordinates  and 
generalized  speeds  must  be  found.  That  is. 


0 

C(J) 

0 

00  = 

c 

0 

S <j) 

- 

9 

1 

0 

C0 

y 

= [u]p 


(2.14) 


where  [u]  is  the  matrix  transformation  from  the  Euler  angular  rates  to  the  absolute 


angular  velocity  vector  for  the  platform,  CO(  . The  derivation  of  matrix  [u]  is  provided  in 
Appendix  A. 

When  substituting  equation  (2.14)  into  (2.13),  the  equation  yields  the  velocity  of 
point  Pj  in  the  { W}  coordinate  system  as  a function  of  the  generalized  coordinates  and 
generalized  speeds.  The  component  of  this  velocity  in  the  (sj  direction  is  the  strut  i 

extension  rate,  which  can  be  obtained  by  dotting  the  velocity  of  point  P,  with 
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the (sj) .vector  as  shown  in  equation  (2.15).  This  operation  results  in  the  magnitude  and 
sense  of  the  strut  extension  rate. 

4 =(vp  - »3  )f  (2.15) 

In  addition  to  the  equation  (2.10),  the  location  of  point  Pj  can  also  be  found  as 
(wr,)f  =(LS,+-r4)f  (2.16) 

Differentiating  equation  (2.16)  yields  the  velocity  of  point  Pi. 

(\  p ) . = (is3  + w , x Ls  3 ),  (2.17) 


where  the  total  strut  angular  velocity  in  {W}  coordinates  is  specified  as 

= K + K + K„  KsJ, 


(2.18) 


Since  (s3  xs3).  =0,  (\  )t  can  also  be  written  as 

(y  „)>  - ^s3  + cos  x Ls3  (2.19) 

and  the  strut  angular  velocity  component  (o  s).  without  the  angular  velocity  about 

A 

S3  axis  is 


(«,),-  = S2  + 0)Sin  s23 )_ 

which  can  also  be  expressed  as 


It  should  be  noted 


(2.20) 


(2.21) 


CQS.  -s3  = 0 

Hence,  a physical  setup,  the  struts  are  typically  designed  to  restrict  free  spin  rotation 
about  its  longitudinal  axis  (s3 ). , therefore,  the  dot  product  of  (c5 s. ). and  (s. ) are  zero 
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( (c5  v • s3 ).  =0).  Thus,  the  strut  angular  velocity  about  (s3 ).  axis  is  (to  v ) and  the  strut 
entities  Body  Bj  and  Body  Ei  have  the  same  angular  velocities  (oy  ) =(00 ; ) = (d>v ) • In 
actuality,  depending  on  joint  designs  for  joint  bj  and  Pj  there  is  small  rotations  about  (s3X 
but  the  dynamics  due  to  Sts  in  the  s3  are  considered  negligible.  Hence,  from  this  point 


forward,  co E and  co/(  are  assumed  equal  to  to  v . 

The  notation  for  the  absolute  locations  and  velocities  of  center  of  mass  of  Bodies  E, 
and  Bi  and  the  platform  joint  Pi  can  be  expressed  as: 


(»rE),  =(Z£s3  + “rb), 

(2.22) 

(WrB),  =(^s3  + “e), 

(2.23) 

(V  J,  = (4*3  + x (4 

A)) 

(2.24) 

(vB),  =Kx(laS3)),. 

(2.25) 

(VJ,.  =vc+(«c.xwrpc) 

, =(AS3+©jx(£Sj))j 

(2.26) 

Strut  angular  accelerations  are  determined  from  the  linear  velocities  and 
accelerations.  To  derive  the  equations  of  motion  in  terms  of  only  the  platform’s  motion, 
the  velocities  and  accelerations  of  the  length  of  the  strut  variable  (L)  can  be  defined  as 
follows: 


From  equation  (2.15), 

A = A;,  = (v  r 'O, 


(2.27) 


L = Z = — (v  • s ) = — (v  ) • (s  ) + (v  )•—(§) 

1 dr  p dr  v,i  v v rJi  dr  ' 


d 


= ^(v‘+((0tXr,,‘)')^’)'  +(VJ,  ■(“  xO, 


(2.28) 


(9Z'Z) 


\s'vrJZ-( 


"Bx£f  T 


= '('»)= r») ='(’“)• 


(sex) 


'(  »7  + 'o)jz) 


((( ' ?7) x 'to)x  '<o+(eS7)x  '»+(!sx  ''<0)73  + £s j)x  '(£s)  = '^bx  ls) 

£'3‘I  ‘ '(CS)  qjIAV  | ^b)  SUTSSOJO 
Xq  paAusp  sq  ubo  puB  sixibs  sqj  sjb  saipoq  ssaqj  joj  suoubj3|30ob  ,iB|n§xiB  aqj  Pxxibs  aqj 
3jb  g pxre  3 ssipoa  aqi  jo  ssppopA  jrujs  JBjnSxiB  Suipuodssxioo  aqj  pup  qoB3  joj  aoxiig 

(p£’Z)  (( £?  "7  x "to)  x "to  + £s'"7  x "»)  = '( "b) 

(( W x Vf0) x w ' S 7 x » + ( £s  x 'to)  '3rji  + ‘s’ '7)  = 

^ ^ V7  x "to  + V7)  — = '(  "b) 

((( l?7)  x ‘ to)  x ‘ to  + ( fS7)  x ' » + ( £s  x 1 0)73  + £S7)  = 

ip 

|£S7x  to  + CS7^  — = 

( ( ,J  .„)  X co)  X co  + | ’"j  ^ X ’»)  + b = '(  "b) 

‘'3'!  '(9Z'Z)  °1  (PZ'Z)  sxioipmba  SutjBqusjsjjip  Ajdxxns  qjiM  paxnxxxrapp  aq  xibo  uuojnqd 
SuiAoiu  aqj  pxre  a puB  3 ssipoa  jo  ssbxxi  jo  isprao  aqj  jo  xioipuspooB  JBauq  aq3 


(Z£'Z) 


(l£'Z) 


'(?>)  + '(‘B).[('(-4.)x  *®)x  '(D+H-J.jx  ’»)+  -’B]=  ^7=  '7  V 

y k'v 

l-s— 


!/  (Zl"  z‘.v  /\\  /II  V II  I/S  \ 

CO+  _ CO/  1 = I <d|  = ( 9J9qM. 


(0£7)  X'®)=  3'Hh  '(*  ’,®+I  ‘'®)  = '((£sx '*»)■( 'gx '<*>)) •'• 

(63‘3)  (( -) x “*'®  + “g x &'®)  = J E? xf  £Zs “‘‘to  + zs "‘to)]  =)£sx  ‘to)  /. 

'(( £S  x ®)  •( £s  X to))  7 + '( £s)  •[(  ( „)  X ’to)  X to  +( '(  ’"j„)  X ’»)  + ’b]  = 

'(£sx  ‘to)-  |(lS7  x to)  + £S7)+  '([s).[(  (’'J,)  x ’to)  x ’to +Q  ’,/j  „)  x ’»)  + ’b|  = 
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The  angular  acceleration  of  the  moving  platform  can  be  obtained  by  differentiating 
to  given  in  equation  (2.14). 

a ,=[U]p  + [U]p  (2.37) 


where  [U]  = 


0 -s  </></>  (c</js0<?)  + S0C00) 

0 c(j)(j)  (sifts  6(f)  - c(f>cOd) 

o o -see 


(2.38) 


2.2.3  Kane’s  Equation 

Using  the  generalized  coordinates  defined  in  equation  (2.1),  Kane’s  method  is  now 
used  to  derive  the  equations  of  motion.  Its  generic  form  is: 


f dv , ( ^ do  ■ 


*=-,  uy  j ( *)  " oca,  ( 

z +z 

M V °Qk  j 7=1  V oqk 


= 0 


k=l,2,...,n 


(2.39) 


d\  f 8(0 . 

where  the  — — and  — - are  the  partial  derivatives  of  the  body  j’s  absolute  linear  and 
dclk  5qk 

angular  velocities  with  respect  to  the  nth  generalized  speed  qn , respectively.  Here  w is 

the  number  of  rigid  bodies  and  n is  the  number  of  the  generalized  coordinates.  The  Fjis 
the  resultant  external  force  acting  on  each  body  j,  while  Fj*  is  the  inertial  force  of  body  j. 
Similarly,  Tj  and  Tj  represent  the  resultant  external  torque/moments  and  inertial  moment 
of  body  j about  its  center  of  mass,  respectively.  Due  to  equation  (2.39)  being  summations 
of  vector  dot  products  and  there  are  n generalized  coordinates,  the  Kane’s  method  results 
in  a set  of  n scalar  equations. 

For  the  Stewart-Gough  PKM,  there  exists  one  moving  platform  (plat),  six- 


extendable  strut  components  (Body  E),  and  six-base  strut  components  (Body  B).  Thus, 
rewriting  Kane’s  equation  for  the  Stewart-Gough  PKM: 
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I 


dqk 


- R - f;  } 


+s 


0(0 


dqk 


— -{t,  -t;} 


+Z 


0V 


dq 


+E 


0(0 


5q 


tMt.-t;} 


°qt  dq„ 


k=l,2,...,6 


(2.40) 


where  (v£,co£)  and  are  the  absolute  (linear  velocity,  angular  velocity)  of 

Body  E and  Body  B,  respectively. 

The  resultant  external  forces  for  each  body  type  are  listed  below. 

F = mpg  + Fw  (2.41) 

F£=wg  + (F  -Ff-FJs3  (2.42) 

F«  = (^V  + F*ib  ~ Fa )s3  + (2.43) 


where  mp  is  the  mass  of  the  platform  and  g is  the  gravitational  field  vector.  Fw  is  the 

applied  tool  and/or  payload  force,  such  as  the  reaction  forces  generated  by  deburring  or 
other  cutting  operations.  The  actuator  force,  Fa,  acts  along  the  vector  s3  with  a 

magnitude  and  sense  of  Fa,  however,  the  friction  force  (F;s3)  and  the  ball  inertial  screw 
forces  (Fm/bS3)  act  in  the  opposite  direction.  The  Fm/b  are  part  of  the  original  unmodeled 
set  of  dynamic  terms  but  are  included  here  for  enlightment.  For  the  external  force  acting 
on  Body  B (FB),  the  forces,  (F^s3,  Fm/bS3,  and  Fas3),  have  the  same  magnitude  but  in 

equal  and  opposite  directions  when  compared  to  those  acting  on  Body  E ( F/r ).  For  all 
bodies,  gravity  forces  are  defined  by  m g terms.  In  the  simulation  so  far,  friction  and 
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ball  screw  forces  are  assumed  to  be  zero.  Since  the  friction  forces  in  the  spherical  and 
Hooke  joint  are  assumed  to  be  zero  and  the  motor  torques  are  modeled  as  equivalent 
linear  actuator  forces,  there  are  no  other  applied  external  torque/moments  to  the  system. 
I.e., 

T =T£=T>0  (2.44) 

For  the  inertial  terms  ( F* , T* ),  it  can  be  written  as: 

F/  = mj*j  (2.45) 

T/  =Oy°Ij+0>_/x(I'offlJ)  (2.46) 

where  a ; is  the  absolute  linear  acceleration  of  the  center  of  mass  of  body  j;  to  . and  a . 

are  the  absolute  angular  velocity  and  acceleration  of  body  j,  respectively;  I"  is  the 

dyadic  term  for  the  inertia  of  the  jth  body  about  its  center  of  mass;  ‘ ° ’ denotes  dyadic 
product  of  the  dyadic  term  with  a vector;  and  ‘ x ’ denotes  the  cross  product  of  vectors. 
The  reader  is  referred  to  reference  [40]  for  detailed  description  of  dyadics. 

The  partial  derivative  terms  from  Kane’s  equation  for  the  Stewart-Gough  PKMs 
are  defined  next.  Using  vector-matrix  notation  and  partition  matrices,  these  derivatives 
can  be  written  as  follows. 


5v 


= [!]■ 


5v 


± = [0]> 


5co 


- = [»]' 


5co 


i = [U] 


(2.47) 


dxc  5p  " J 5vc  L J 5p 

where  [I]  is  a (3x3)  identity  matrix  and  [0]  is  a (3x3)  zero  matrix.  For  each  strut  i,  the 


following  equations  are  stated  with  subscript  i omitted  for  brevity.  Note,  i , j,  k are  the 
unit  vectors  of  {W},  Xw,  Yw,  Zw  respectively.  The  ; in[  : ] are  used  in  forming 


partitions  in  the  [ ] matrix. 
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(2.52) 


d«*s 


d(oR  do 


ap  ap 


ap^  = (l )fcxtuW>x’r^))  ; Mu]c< 


Col  .2 


J i 


(2.53) 

where  [U]Co/  / is  the  ith  column  of  matrix  [U], 

2.2.4  Equations  of  Motion  (Force-Based  Model) 

Once  the  above  equations  are  substituted  into  the  Kane’s  equation,  a set  of  six 
scalar  equations  will  result,  where  the  six  (FJj  terms  can  be  factored  out.  The  force- 

based  model  uses  the  resulting  model  expression  for  the  vector  of  these  six  actuator  force 
terms  ( Ffl  )j.  When  expressing  these  six  equations  of  motion  in  an  ordered,  compact 

vector-matrix  form  per  the  defined  generalized  coordinates,  the  coefficients  of  these  (F„)i 
terms  yield  a pseudo-Jacobian  matrix 


30 


J = 


(2.54) 


analogous  to  the  traditional  Jacobian  definition  for  PKMs,  refer  to  Appendix  B (equation 
(B.33)-(B.36))  for  further  derivation  details.  That  is,  the  Pseudo-Jacobian  matrix  can  be 
physically  interpreted  as  the  pose  of  machine  expressed  in  the  Euler  angles,  especially  if 
one  takes  into  account  that  the  [u]  matrix,  defined  in  equation  (2.14),  is  the 
transformation  matrix  from  the  Euler  angular  rates  to  the  absolute  angular  velocity  of  the 
platform  in  { W}  coordinates.  Solving  for  the  actuator  forces  yields. 


F 
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= J ' 
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L ext  J 

T 

L m/h  J 
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(2.55) 


where  F(,  denotes  a vector  of  the  magnitude  and  sense  of  the  six  actuator  forces,  each 


applied  along  their  respective  strut  axes  (s3  )i.  In  the  above  equation  of  motion,  equation 

(2.55),  the  dynamic  forces  and  moment  contributions  and  the  external  loads  are  grouped 
by  their  source.  Details  of  the  components  of  these  forces  and  moments  are  expressed  in 
Appendix  B.  For  the  platform  acceleration  dependent  terms,  the  forces  and  moments  are 

grouped  as  [f^  T;m  ] . The  gravitational  dependent  terms  for  the  platform  are 

[ y 7'  IT’ 

TKrp  J . Similarly,  strut  acceleration  and  the  gravitation  dependent 

terms  in  equation  (2.55)  are  grouped  in  Fvn/  T,/]  and  [f^/  T / T, 
respectively.  The  Coriolis  and  centripetal  acceleration- like  terms  are  grouped  in 

7'  ^ ^ 7'  |7’  I j j T 

Fw  Tcor  J and  [Fcai  Tcen  , respectively,  where  they  contain  terms  that  are  a 
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function  of  strut  and  platform  velocities  products.  As  expected,  the  resulting  equations  of 
motion  are  highly  nonlinear.  Other  sources  of  nonlinearities  are  joint  frictions  which  are 

[ J1  y |7' 

F^.,  T I and  are  dependent  on  the  PKM  instantaneous 


configuration,  velocities  and  friction  modeled.  The  [f^/  T frJ  ] terms  are  set  to 
zero  for  the  test  without  friction,  where  the  friction  modeling  (Chapter  4)  is  added  for  the 

[ Y y Y 

F ^ T ,xl  J consists  of  the  external  loads  acting  on  the  PKM. 

[ Y y P 

F„,/A  T m/h  J represents  the  terms  due  to  the  motor  and  ball  screw  dynamics. 

The  [¥umd  Tumd  ] denotes  the  collection  of  the  other  unmodeled  dynamics  term.  In 


the  remainder  of  this  dissertation,  the  [Fm//  Tm/J  ] and  [F;h;/  T J J terms  are 

assumed  to  be  negligible.  Table  2. 1 shows  a summary  of  the  components  associated  with 
each  term  of  equation  (2.55).  For  further  details  on  generating  the  equations  of  motion, 
the  reader  is  referred  to  Walker  [39]  and  Register  [36]. 

2.2.5  The  Energy  Method 

The  prior  sections  in  this  dissertation  were  dedicated  to  obtaining  the  full  set  of 
dynamic  equations  of  motion  to  formulate  the  model  basis.  For  the  parameter 
identification  algorithms,  the  implementation  of  the  identification  algorithm  with  the 
dynamic  equations  of  motion  as  the  model  basis  involves  measuring  position,  velocity, 
acceleration  and  actuator  forces  as  the  inputs  to  the  identification  process.  Measurements 
of  velocity  and  acceleration  are  typically  prone  to  errors.  When  data  is  not  directly 
measurable  but  must  be  calculated  numerically  from  other  measured  data,  e.g.,  position 
data,  errors  can  rapidly  grow  in  magnitude  especially  with  each  level  of  differentiation 
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that  is  calculated  using  data  containing  noise.  To  circumvent  this  problem,  one  can 
instrument  the  setup  with  better  quality  sensors,  which  tends  to  drive  up  costs  and  require 
greater  engineering  expertise.  An  alternative  would  be  to  use  as  the  model  base  the 
energy  method  suggested  by  Bhattacharya  [33],  In  doing  so,  one  does  not  need  to 
measure  accelerations,  thus  eliminating  one  level  of  noise  amplification  in  calculated 
inputs  to  the  identification  algorithm.  Thus,  an  energy-based  model  approach  is 
investigated  as  a potential  alternative  to  the  equation  of  motion  based  approach,  referred 
to  as  the  force-based  approach  in  this  dissertation. 

For  the  energy-based  method,  the  PKM  is  assumed  to  move  in  a singular-free 

trajectory,  permitting  an  energy  difference  determination  for  time  step  through 

(X/t  + l’Xjt  + l ) . The  energy  equations  developed  by  Bhattacharya  [33]  are  defined  as  the 
following  scalar  function. 

(x*+. > > w*+, ) - K (**  , ** , wk ) = f F • x dt  = £ (F  ) (AZ,  )k  (2.56) 

i=l  ' 

where  hk  (xk  ,xk,wk)is  the  total  energy  (kinetic  plus  potential  energy)  at  the  kth  states, 
K+i(xk+\’Xk+\’wk+i)  is  total  energy  (kinetic  plus  potential  energy)  at  the  (k+l)th  state,  wk 
and  wk+x  are  the  current  kth  and  (k+l)th  inertia  parameter  estimates.  [Fa  ) is  the  ith  strut 
force  at  the  kth  state,  (AL  is  the  change  in  the  ith  strut  length  at  k'h  state,  n is  the  number 
of  struts. 

For  the  PKMs  modeled  in  this  dissertation,  the  mass  matrix  is  the  same  as  that 
defined  in  Appendix  B (equations  (B.2)  - (B.8)  and  (B.25)).  The  energy  at  kth  state  is 
expressed  as  follows. 


hk(xk,xk,wk)  = Tk+Vk 


(2.57) 
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where  the  kinetic  energy,  j 


is  a scalar  quantity  on  dependent  of  the 


choice  of  generalized  coordinates  as  illustrated  in  Appendix  B,  equations  (B.26)-(B.32) 
and  the  total  potential  energy  in  the  system  is  given  by 


K = (-1) 


T W 

mg  r 

p O c 


+ Z['»  g ' “ r, 


l+Zkg"r„ 


L /-I  i-1 

2.3  Parameter  Identification  Algorithms 

In  this  section,  the  proposed  methodology  for  identifying  the  dynamic  system 
parameters  will  be  shown.  In  general,  the  minimal  set  of  dynamic  parameters  consists  of 
ten  parameters  per  rigid  body:  the  mass,  inertia  tensor  and  the  locations  of  the  center  of 
mass  per  rigid  body  in  the  system.  These  terms  can  be  extracted  from  the  system’s 
theoretical  model-base  equations  forming  the  set  of  potentially  identifiable  parameters. 

As  the  system  equations  become  more  complex  and  nonlinear,  the  extraction  of  the 
parameters  in  forming  the  traditional  vector  set  of  identifiable  parameters,  as  required  by 
least  squares  (LS)  and  extended  Kalman  filter  (EKF)  methods,  becomes  insurmountable. 
Since  the  ultimate  goal  for  this  dissertation  is  to  develop  a methodology  for  determining 
accurate  dynamic  system  parameter  values  for  parallel  kinematic  machines  (a  highly 
nonlinear  system),  the  ease  of  implementing  the  estimation  process  is  critical.  The 
system  identification  process  for  parameter  estimation  will  be  chosen  between  the 
following  methods: 

• Least  square  (LS) 

• Extended  Kalman  filter  (EKF) 

• Unscented  Kalman  filter  (UKF) 
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2.3.1  Least  Squares  Method  (LS) 

Since  the  parameter  estimation  algorithm  has  a long  history,  the  pros  and  cons  of 
the  LS  method  are  well  known  to  the  readers.  For  linear  systems,  many  researchers 
[23,32]  have  used  the  LS  method  for  system  identification.  The  system  identification 
process  with  a linear  trajectory  and  linear  parametric  modeling  is  successful  for  most 
cases.  However,  when  nonlinear  terms  are  present,  the  nonlinearity  of  the  system  may 
cause  defects  in  the  identification  analysis  and  lead  to  convergence  problems. 

The  LS  method  is  based  on  the  equation  (on  kth  iteration  step) 


where  reference  residual,  e = (yk  - dt ) . , the  yk  is  the  updated  output  (calculated  from 
measured  states  (xk)  and  estimated  parameters  (wk ))  for  kth  iteration  step,  the  d,  is  the 

measured  (desired)  output,  and  m is  the  number  of  residual  equations. 

The  linearized  form  yields 


where  wko  is  the  nominal  value  for  parameters  at  kth  iteration  step,  e I is  the  ith 

' 1 wko  j 

residual  evaluated  at  Wko , de,  , is  the  partial  derivative  of  the  jth  parameter  and  is 


number  of  parameters.  For  m equations  or  poses,  the  vector  equation  form  is 


Min / = £>,. 2 


(2.58) 


(2.59) 


dw 


evaluated  at  wko , the  is  the  change  in  the  nth  parameter  for  j=l,n,  and  n is  the 


e*  = k +J,.s,k^k  + H.O.T. 


(2.60) 
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J LS,k  ~ 


dex 

dw. 


de„ 


de] 

dwn 

de... 


dw , dw.. 


(2.61) 


and  H.O.T.  stands  for  higher  order  terms.  If  ignoring  the  H.O.T. , the  objective  function 


is 


Minf  = (ek  |no  +Jls,kAwk)T(ek  \ +JLSJ[Awk) 


(2.62) 


The  minimum  value  is  derived  by  setting  - ^ = 0 together  with  ^ > 0 . 


dAw. 


dAw,, 


JL 

dAw, 


2-^ LS,kek  lwto  +(JkS,k J IS, k + ^LS,k^LS,k)^Wk  ~ 0 


(2.63) 


which  yields 


Awk  (J l.S,kJ LS,k  ) J LS,ke 


* lw*„ 


(2.64) 


and 


wk+ 1 =wk+ 


(2.65) 


The  LS  program  iterates  until  the  largest  value  in  A wk  is  less  than  desired  tolerance  to 


obtain  the  estimation  of  W,  . 

k 

2.3.2  Extended  Kalman  Filter  (EKF)  and  Unscented  Kalman  Filter  (UKF) 

The  EKF  method  has  widely  been  used  because  of  its  benefits  [36],  For  the  PKM 
system,  the  EKF’s  ability  to  handle  noisy  data  and  highly  non-linear,  multi-input/multi- 
output systems  is  attractive.  However  recently,  Wan  and  Van  der  Merwe  [38] 
demonstrated  the  UKF  as  an  alternative  to  EKF,  pointing  out  the  flaws  and  shortcomings 
of  the  EKF.  The  EKF  involves  a linearization  approximation  within  its  algorithm  when 
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applied  to  nonlinear  systems.  This  causes  an  elimination  of  the  expectations  for  the 
optimal  prediction  terms,  which  induces  large  errors  in  the  true  posterior  mean  and 
covariance.  Thus,  Wan  and  Van  der  Merwe  [38]  proposed  the  UKF  method  to  overcome 
these  algorithm  limitations.  Instead  of  using  a first-order  linearization  approximation,  the 
UKF  uses  an  optimal  prediction  of  the  input.  The  technique  that  eliminates  the  first-order 
approximations  is  called  Unscented  Transformation  (UT).  The  UKF  estimates  in  two 
steps.  The  first  step  evaluates  and  estimates  the  existing  value  of  the  initial  guess  of  the 
parameters  using  the  unscented  transformation.  In  the  second  step,  the  parameterized 
vector  is  estimated  using  a Kalman  gain  (weighting)  technique.  The  following  section 
details  the  theory  of  the  EKF  and  UKF  and  further  highlights  the  main  differences 
between  EKF  and  the  UKF. 

2.3.2. 1 Extended  Kalman  Filter  Parameter  Estimation 

For  the  basic  state-space  estimation,  recursive  estimation  for  input  xk  can  be 
expressed  as 

xk  =(prediction  of  xk  )+&*•[  yk  -(prediction  of  yk )]  (2.66) 

where  kk  is  the  optimal  gain  term.  The  optimal  terms  for  this  recursion  are 

**  =TO-i.Ui)] 

k*=P.,„K„  (2.67) 

h = £[«(*»,»»)] 

where  the  optimal  prediction  of  xk  is  expressed  as  xk , and  the  optimal  prediction  of  the 
yk  is  expressed  as  yk  . The  vkA  represents  the  process  noise,  nk  stands  for  the 
measurement  noise.  The  E[f(  )]  denotes  the  expected  value  of  the  function  F(  ) and  the 
F[H(  )]  denotes  the  expected  value  of  the  function  H(  ).  The  determination  of  this 
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expectation  function  depends  on  the  estimation  procedure  selected.  The  P function  is 

xk)>k 


defined  as  the  covariance  of  the  x-  and  yk  functions.  The  Kalman  filter  and  extended 

Kalman  filter  techniques  approximate  these  quantities  using  linearization,  i.e., 

**  ~ F(xk_„v) 


* Vw 
h » Wirt 


(2.68) 


These  approximations,  however,  can  lead  to  large  errors  [38]  when  applied  to  nonlinear 
systems. 

2.3.2.2  UKF  Parameter  Estimation 

The  UKF  finds  its  roots  in  the  1995  development  of  the  Unscented  Transformation 
(UT)  used  for  filtering  nonlinear  systems  [42,  43],  Julier  and  Uhlmann  [44]  later  applied 
the  UT  to  extend  Kalman  filtering  in  the  development  of  a generalized  method  for 
estimation  of  nonlinear  systems.  The  classic  nonlinear  mapping  takes  the  form. 

yk=G(xk,wk)  (2.69) 

where  xk  is  the  measured  states  and  yk  is  the  updated  output  (calculated  from  measured 
states  (xk)  and  estimated  parameters  ( wt )),  and  the  nonlinear  mapping  G is  parameterized 
by  wk . The  error  of  the  mapping  is  ek  =dk  - G(xk , wk ) , and  the  goal  is  to  minimize  the 
expected  squared  error.  The  dk  stands  for  the  measured  (desired)  output.  The  estimation 
procedure  for  the  parameter  can  be  defined  by  the  following  expressions. 

(2  70) 

dk  =G(xk,wt)  + e. 


where  wk  corresponds  to  a stationary  process  with  an  identity  state  transition  matrix, 
driven  by  the  noise  uk . The  output  yk  corresponds  to  the  nonlinear  observation  on  wk . 
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The  wk  parameters  are  updated  via  equation  (2.65)  for  each  kth  iteration.  The 

distinguishing  difference  of  UKF  from  other  methods  is  that  it  uses  sigma  points  and  it  is 
efficient,  derivative-free  for  on-line  learning. 

23.2.3  UKF  Parameter  Estimation  Procedure 

The  UKF  addresses  the  problems  faced  in  using  of  the  EKF  when  applied  to 
nonlinear  systems.  Using  a carefully  chosen  set  of  sigma  points,  the  UKF  captures  the 
true  mean  and  covariance  of  the  Gaussian  Random  Variables  (GRV).  A method  for 
calculating  the  statistics  of  a random  variable  is  the  unscented  transformation  (UT).  The 
UKF  is  a direct  extension  of  the  UT.  First,  one  assumes  the  parameter  ( wk ) to  be 

estimated  has  mean  wk  and  covariance  P . To  calculate  the  statistics  of  the  output  yk 


in  UKF,  takes  the  form  yi  k/k_x , a matrix  Wjklk_x  (with  corresponding  weights  W.0'0)  is 
formed,  according  to  the  following: 


= wk  i = 0 

= *.+(,/(!  +A)P.,  ),=*,+ 

(yJp7)i-+i  = \...,L 

= *,-W(Z +w.  )»=*.- 

w ("') 

= A/(Z+A) 

iV’ 

= A /(Z,  + A)  + (1  — CC~  + y9) 

(2.71) 


(2.72) 


W.(m)=w,(c)  = 1 /{2(L  + A)}  -»  z = l,...,2L 

where  A = or  (L  + k)  - L is  a scaling  factor.  The  a determines  the  spread  of  the  sigma 
points  (/F  t/i  i ) around  Wk  and  is  usually  set  to  a small  positive  value  (e.g.  0.001).  The 

L stands  for  the  dimension  of  the  parameter.  The  k is  an  additional  scaling  factor  and  is 
usually  set  to  zero.  The  /?  is  used  to  incorporate  prior  knowledge  of  the  distribution  of  x 
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(for  Gaussian  distributions,  fi  =2  is  optimal).  The  (yj(L  + A)PW ) . is  the  ith  row  of  the 
matrix  square  root.  The  sigma  points  ( W.  k /k  l ) calculation  is  the  checkpoint  of  whether  the 
current  parameter  values  are  in  the  convergence  range  or  not.  This  calculation  is  shown 
on  Figures  2.8  and  2.9  via  block  diagram  of  the  UT  illustrating  the  generation  of  d k and 


Pj  u ’ an^  a fl°wchart  of  the  UKF  algorithm.  The  mean  and  covariance  for  output  d , are 
approximated  using  the  weighting  sample. 


4=2X"V, 


k/k-l 


(2.73) 


p-j-di  =EW/  -dk][yiklk_,  -dkY  +r; 

/'= 0 


(2.74) 


In  UKF  procedure,  there  are  several  important  points  to  check  to  use  this  method 
[45],  First,  the  convergence  rate  and  tracking  performance  are  directly  influenced  by  the 


choice  of  the  covariance  E 


T \ r • r 

ukuk  J=  Rk  . Units  for  Rk  are  the  same 


as 


WkWK 


Basically,  as  Rkr  gets  larger,  the  parameters  are  converging  quickly.  To  choose  Rkr 

value,  there  are  several  options.  One  option  is  the  randomly  choosing  a diagonal  matrix. 
Another  option  is  the  Robbins-Monro  stochastic  approximation,  which  is 

Rkr  ={\-a)R[_x  +aKk(dk -G(xk,w))(dk -G(xk,w))T  (K?)T  (2.75) 


where  Rkr  is  also  restricted  to  be  a diagonal  matrix.  Wan  et  al  [45]  also  found  that  a 
Robbins-Monro  stochastic  approximation  method  converges  better,  in  general.  Thus,  in 
this  dissertation,  equation  (2.75)  was  employed  to  calculate  and  update  Rkr . Secondly, 
complete  convergence  can  be  achievable  when  the  covariance  PWt  goes  to  zero. 
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Determining  Pw*  will  produce  a good  basis  for  excellent  convergence.  In  Chapter  3,  the 

nature  of  the  sigma  points  and  demonstration  of  UKF  convergence  is  further  explained 
using  the  example  of  an  inverted  double  pendulum. 

2.4  Implementation  Issues 

From  prior  research,  the  base  and  moving  platform  joint  locations  and  joint  location 
errors  were  found  to  have  an  impact  on  the  overall  system  dynamics  (Johnston  [35], 
Walker  [39]).  Since  the  main  goal  of  this  dissertation  is  estimating  the  dynamic 
parameters  of  the  system  understanding  the  different  contributions  of  the  different 
dynamic  groups  defined  in  equation  (2.55)  would  provide  insight  as  to  which  parameters 
may  have  a higher  probability  of  being  identifiable  as  well  as  those  that  may  not  be 
identifiable.  This  section  presents  a study  of  the  above  PKM  dynamics  grouped 
according  to  inertial  and  gravitational  components  for  three  different  Stewart-Gough 
PKM  configurations  by  comparing  the  relative  magnitude  of  each  term,  an  observation  of 
the  relative  contribution  to  the  dynamics  response  can  be  made.  Based  on  this 
observation,  a prediction  of  the  dominant  parameters  can  be  made.  In  addition,  this  study 
provides  insight  as  to  the  existence  of  trajectories  that  will  yield  sufficiently  dynamically 
rich  data  necessary  for  parameter  identification.  The  specific  platforms  included  in  the 
study  are  the  University  of  Florida  Special  6-6  PKM,  the  Hexel  6-3  PKM  machine  tool, 
and  the  Ford  General  6-6  PKM  vehicle  simulator,  shown  in  Figure  2.4  and  described  in 
section  2.1.1.  Their  analysis  also  provides  insight  to  how  geometry  (e.g.,  joint  locations) 
plays  a role  in  the  dynamics.  Comparative  observations  are  made  using  Euclidean  norms 
and  sample  trajectory-based  characterizations. 
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2.4.1  The  Characterization  of  PKM  Dynamics  in  Terms  of  the  System  Identification 

2.4.1. 1 System  Model 

In  order  to  eliminate  the  overall  PKM  size  effect  on  the  results  and  to  make  a more 
direct  comparison  between  the  PKM  variants  presented  in  section  2.1.1  and  Figure  2.4, 
their  relative  sizes  were  scaled  to  a physical  size  of  the  same  order  of  dimensional  scale 
as  the  UF-PKM.  That  is,  the  moving  platforms  of  each  PKM  were  scaled  to  form  a 
theoretical  equivalent  equilateral  triangle  the  same  size  as  UF-PKM’s  moving  platform. 
The  Ford  equivalent  resulted  in  truncated  corners  of  the  equivalent  triangle.  All  other 
dimensions  maintained  the  same  relative  geometric  spacing  given  by  side: SIDE.  The 
size  of  the  UF  PKM  is  SIDE  = 0.8128  m and  ratio  side:SIDE  = 0.750.  The  Hexel  PKM 
has  the  joint  configuration  of  the  PKM  shown  in  Figure  2.4  (b),  size  scaled  to  SIDE  = 
0.4722  m and  ratio  side:SIDE  = 0.6096.  Physically,  it  is  constructed  with  platform  joints 
paired  in  bifurcated  ball  joints  where  each  strut  is  the  screw  of  the  ballscrew  rigidly 
attached  to  the  bifurcated  ball  joint  at  one  end,  and  passing  through  the  center  of  the  base 
vertex  (motor  driven  ball  nut).  The  Ford  PKM  has  a near  octahedral  layout  similar  to  the 
Hexel,  but  with  the  platform  joints  in  close  pairs,  as  the  base  joints.  The  Ford  PKM  has 
the  joint  configuration  of  the  PKM  in  Figure  2.4  (c),  size  scaled  to  SIDE  = 0.5276  m and 
ratio  side:SIDE  = 1.12. 

2.4. 1.2  Trajectory  Selection 

In  order  to  capture  the  dynamics,  a motion  must  be  specified.  Unfortunately,  the 
resulting  dynamics  is  very  dependent  on  the  choice  of  motion.  Hence,  for  the  study 
presented  here,  various  trajectories  were  selected.  Sinusoidal  translation  trajectories  of 
the  platform  in  each  of  the  X,  Y,  Z and  X-Y  directions  about  the  ‘home’  position  of  the 
platform  (center  of  its  workspace  at  mid-stroke  of  struts)  were  selected  and  run  for  low, 
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medium  and  high  velocities  of  the  platform.  Similarly,  sinusoidal  trajectories  in 
rotational  variables,  Euler  (<|),  0,  \\i),  of  the  platform  about  its  home  position  were  selected 
and  run  at  3 speed  levels.  For  each  of  the  above  trajectories,  a platform  offset  in  the  X-Y 
plane  from  the  home  position  was  also  introduced  for  additional  evaluation.  The 
amplitudes  of  the  motion  were  selected  so  that  the  PKM  would  remain  within  its  physical 
workspace  but  span  a reasonable  breadth  (50%  or  greater).  Table  2.2  lists  a subset  of 
these  trajectories  by  the  Case  Numbers  used  in  Figures  2.10  through  2.15. 

The  Hexel  and  Ford  PKMs  were  simulated  for  the  same  trajectory  scaled  up  such 
that  each  machine  moved  from  the  center  of  their  respective  workspace  to  within  20%  of 
the  work  boundary.  Thus,  while  the  UF  PKM  moved  to  a final  position  of  X=Y=0.04  m, 
the  Hexel  PKM  moved  to  X=Y=0.3  m and  the  Ford  PKM  moved  to  X=Y=0.6  m. 

Velocity  and  acceleration  profiles  were  also  scaled  according  to  each  PKMs  maximum 
feed  [UF  + 0.3m/s  (4m/  sec2),  Hexel  ^ 0.3  m/s  (3  m/sec2),  Ford  1 m/s  (10  m/sec2)]. 
Tool  orientations  were  15  degrees  for  Hexel,  10  degrees  for  Ford,  and  5 degrees  for  UF 
PKM.. 

2. 4.1.3  Metrics  of  Identifiable  Terms  and  Effort  Expended 

In  order  to  establish  which  dynamics  terms  belong  to  an  identifiable  set.  Euclidean 

norms  of  the  Fa  vector  of  actuator  inputs  and  the  model’s  individual  J1  [f7  T7  ] 

vectors  (from  equation  (2.55))  were  calculated  along  each  of  the  above  sinusoidal  motion 
trajectories  (listed  in  Table  2.2).  By  normalizing  the  Euclidean  norms  with  respect  to  the 
mean  value  of  the  actuator  Fa  norm  along  a given  trajectory,  one  achieves  a ‘relative’ 
percentage  measure  of  the  theoretical  dominance  of  each  dynamics  term  in  a 
corresponding  measured  Fa  signal.  Figures  2. 1 0 through  2. 1 5 provide  a representative 
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set  of  the  following  norm  evaluations  along  each  given  trajectory  that  was  used  for 
quantifying  the  results  of  this  study.  For  Figures  2.10  through  2.13: 
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For  Figures  2.14  and  2.15,  the  peak-to-peak  value  of  the  F norm  along  a given 

trajectory  is  calculated  to  provide  a metric  of  the  relative  variations  in  the  force  where  the 
larger  values  indicate  a greater  likelihood  that  parameters  may  be  identifiable. 

Fa  norm  : {(max  - min)  of  F ' J4  mean  {jFjFa  $*100  (2.82) 


pm  norm  : j (max  - min)  of 


sm  norm  : j (max  - min)  of 


cor  norm  : j (max  - min)  of 


cen  norm : < (max  - min)  of 


grp  norm  : < (max  - min)  of 


1. 
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)T J 1 [f  r T rf 

/ L pm  pm  J 
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grs  norm  : j (max  - min)  of 


t/MC  TJIIJ  1F,/  T„r]'  1 + mean  100  <2'88) 


2. 4.1.4  Dynamic  Force  Contributions 

In  this  section,  the  PKMs  shown  in  Figure  2.4  (a)  are  assumed  to  have  the 
following  mass  parameters  of  the  UF  PKM:  mass  of  platform  mp  = 4.82  kg  with 
Ixx=Iyy=0.0908  kg-m"  and  1^=0.181  kg-m2.  The  strut  masses  are  mb=7.72  kg  with 
Ixxb=Iyyb=0.552  kg-m"  and  1^=0.00751  kg-m2,  and  me=1.79  kg  with  Ixxe=Iyye=0.0183  kg- 
m2  and  Izze=0.00160  kg-m2.  The  ratio  mp/(6*me)  = 0.45.  Using  the  prescribed 
trajectories  of  section  2.1,  first  paragraph,  the  corresponding  contribution  of  the 
individual  dynamic  components  for  each  PKM  configuration  is  presented  and  discussed 
in  section  2.1.1. 

For  the  Hexel  PKM,  masses  and  inertias  were  estimated  using  dimensional  and 
material  data.  The  mass  of  the  platform  was  estimated  to  be  93.0  kg  with  Ixx=Iyy=l . 1 8 
kg- nr  and  Izz=1.37  kg-m".  Strut  masses  were  calculated  as  mb=2.1 1 kg  with  IXXb=Iyyb= 
Izzb=0.00135  kg-m2,  and  me=4.64  kg  with  Ixxe=Iyye=  0.556  kg-m2  and  lw  = 0.000362  kg- 
rrf . The  Flexel  PKM  ratio  mp/(6*me)  = 3.34. 

Most  Ford  PKM  parameters  were  obtained  from  measurements  or  design 
specifications.  In  some  cases,  parameters  had  to  be  extrapolated  from  the  given  data. 

The  values  used  in  this  study  are  as  follows.  Platform  mass  was  7500  kg  with  Ixx= 
Iyy=6.42xl03  kg-m“  and  1^=1 .28x1 04  kg-m2.  Strut  masses  were  nib=100  kg  with 
I\xb=Iyyb=85.6  kg-m"  and  IZZb=0.500  kg-m2.  and  me=54.5  kg  with  Ixxe=Iyye=l  .83  kg-m2  and 
lzze=0. 174  kg-m".  The  Ford  PKM  ratio  mp/(6*me)  = 22.9. 
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2.4. 1.5  Identifiable  metrics  parameters 

Implementing  the  series  of  Euclidean-based  norms,  a Matlab  simulation  of  the 
PKM  dynamics  was  run  for  each  of  the  sinusoidal  motions  trajectories.  The  F„  mean 

(stands  for  mean^Fj  ■ Fa ))  values  were  found  to  be  relatively  large  for  each  motion 

case  studied.  Generally,  these  large  means  were  attributed  to  the  gravitational  effects  at 
low  velocities  (0.0254*  n m/sec),  while  at  higher  platform  velocities  were  attributed  to 

the  gravitational  and  the  [F  w 7 T./]  and/or  [F/jb/  T 7 } effects.  Figures  2.10 


and  2. 1 1 illustrate  the  relative  contribution  trend  in  the  dynamics.  The  relative  degree  of 


contribution  between  the 


F,/  Tm7]  and  [F/  T,/  ] shifts  with  the  mass  of  the 


platform-to-mass  of  the  struts  ratio,  mp/(6*me),  as  illustrated  by  Figures  2.12  and  2.13. 


As  this  mass  ratio  increases 


, [f  7'  T 

5 L pm  } 


T 
pm 


l 


relative  dominance  increases.  The  same 


trends  were  found  for  each  PKM  configuration  of  Figure  2.4,  as  well  as  for  the  Hexel  and 
Ford  PKMs.  In  terms  of  geometry  affects,  the  PKMs  of  Figures  2.4  (b)  and  2.4  (c) 
exhibited  the  same  general  trend  shifts  as  those  of  PKM  in  Figure  2.4  (a).  However 
examining  Figures  2.1 1 and  2.13  closely,  one  also  finds  that  as  the  joint  locations 


approach  the  6-3  configuration,  the  amount  of  the  F ' T \ and  F ‘ T | 

L V7w  sm  J L pm  pm  J 

[T  T r 

Fvm  T m J in  the  horizontal  and 

angular  motion  cases. 

Figures  2.14  and  2.15  show  the  results  of  the  force  variation  norms  defined  relative 
to  the  average  of  the  Fa  norm  (equations  (2.82)  through  (2.88))  for  six  of  the  motion 
cases  studied.  As  can  been  seen  from  Figure  2.14,  the  variations  of  all  forces  were 
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relatively  low  for  low  platform  velocities.  Figure  2. 1 5 illustrates  that  not  only  is  the 

[Fm  Tm  ] a dominant  dynamic  term  based  on  the  mean  values  (Figure  2.11),  but 

truly  dominates  in  terms  of  force  variations  relative  to  the  Fa  norm.  The  gravitational 
forces  have  small  variations  in  their  contributions  to  Fa.  For  high  platform  velocities,  the 

T „/  ] exhibits  a modest  contribution  to  the  variations  in  F„.  As  the  relative 


platform  mass  increases,  the  [Fv(/  Tj  } and  [fJ  Tpn'  } switch  roles. 

Although  a very  modest  effect,  Figure  2.14  delineates  that  offset  in  the  motion  generates 
greater  variations  in  both  the  gravitational  forces  and  Fa. 

Contrary  to  the  above  observations,  the  Hexel  PKM  exhibited  relatively  large  Fa 

and  [Fot,  T;  ] variation  norms  (20  to  58)  for  all  motions  studied  where  the  motion 

had  an  offset  from  its  center  Z axis.  This  is  attributed  to  the  geometric  effect  of  the  small 
platform  size  relative  to  the  base  size,  causing  increased  gravity  loading  variations  within 
the  struts. 

To  summarize  and  referring  to  Table  2.1,  the  selected  motions  of  this  study  indicate 
that  the  identifiable  set  of  dynamics  terms  for  extendable  strut  PKMs  will  be 

[F*»  T«,  f ’ k,  T/ ] , [F/  Tj  ] and  [f„,7  T/  J . The  Coriolis 


and  centripetal  dynamics  terms  were  found  to  have  very  small  effects  in  all  cases.  Hence, 
the  dynamic  parameters  within  these  quantities  may  not  be  as  easily  identifiable, 
especially  if  the  F„  signals  are  noisy.  (Note,  for  PKMs  that  have  significantly  greater 
orientation  capabilities  this  may  not  be  the  case.)  In  addition,  changes  in  the  mp/(6*me) 
ratio  and  the  motion  type  and  rates  exhibit  more  apparent  affects  on  the  relative 
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dominance  of  the  dynamic  terms  in  a given  trajectory  than  that  observed  for  the  changes 
in  geometry.  The  above  findings  indicate  that  slow  motions  of  the  platform  can  be 


utilized  to  first  identify  the  masses  and  the  center  of  mass  locations  in  a simplified, 
multiple  step  identification  method.  As  the  platform  size  relative  to  the  base  decreases, 
the  identification  of  the  other  dynamic  terms  however  may  be  more  challenging.  The 
horizontal  motions  (parallel  to  the  base)  are  more  suitable  for  identifying  parameters 


within  the 


F,;  T m ] and  [f^  T ; ] terms.  While  both  vertical  and  slow 


motions  are  dominated  by  the  F ‘ T ‘ and  F ' T ' [ terms 

J L grs  grs  J L grp  grp  J 


2. 4. 1.6  Summary 

To  summarize  the  findings  of  this  section,  an  in-depth  study  on  the  dynamics  has 
conducted  for  three  parallel  kinematic  mechanism  configurations.  From  this  study,  the 
determination  of  which  dynamics  terms  are  ‘candidates’  for  identification  have  been 
established,  highlighting  the  potential  limitations  due  to  similarities  in  individual 
dynamic  contributions  and  other  system  uncertainties.  Consequently,  the  study  indicates 
that  dynamic  system  identification  may  require  an  ‘iterative’  methodology  for  system 

identification,  so  as  to  extract  additional  dynamic  parameters.  In  brief,  when  F^ 
dominates  the  signal,  then  one  can  expect  to  be  able  to  extract  the  value  of  m easily. 
Similarly,  for  when  F.m  dominates  then  one  can  expect  to  be  able  to  extract  m , , mh  ,Ie 


and  lb- 
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Overall,  the  results  show  that  trajectories  sufficiently  rich  in  dynamic  excitation 
(commonly  refer  to  as  dynamically  exciting  trajectories)  are  plausible,  hence  the  ability 
to  identify  the  prominent  dynamic  parameters  from  experimental  data  exists. 

2.4.2  The  Parameter  Sensitivity  in  Sensor  Location 

On-line  parameter  estimation  has  been  shown  to  be  a difficult  task  for  the 
generalized  Stewart  PKM  [33],  Therefore,  in  an  effort  to  decrease  the  complexity  and 
difficulties  of  actually  implementing  parameter  identification  methods,  the  effects  of 
limited  sensing  instrumentation  are  studied  in  this  section.  The  motivation  of  this  study 
as  well  as  the  one  just  prior  to  this  section  was  to  provide  a level  of  understanding  and  to 
illuminate  areas  of  potential  simplification  necessary  for  realizing  an  industrial 
acceptance  of  the  presented  methodology.  Thus,  to  simplify  the  estimation  process,  this 
section  focuses  on  the  validity  of  various  simplifying  assumptions.  That  is,  the  author 
suggests  using  the  geometric  center  in  lieu  of  estimating  the  true  center  of  mass  from  the 
measured  data  for  the  sensor  location. 

2.4.2. 1 Definition  of  Sensitivity 

To  identify  the  dynamic  terms  of  special  6-6  PKM.  it  is  required  to  experimentally 
measure  the  input  torques  and  motion  (e.g.,  motor  current  and  encoder,  respectively). 
Since  the  inaccuracy  of  the  motor  current  causes  significant  noise  and  errors,  new 
measurement  device  techniques  are  sought.  Such  techniques  could  consist  of 
force/torque  sensors  placed  directly  on  the  platform  or  within  the  struts.  For  motion,  the 
use  of  laser  trackers  or  other  optical  devices  are  expensive  and  require  unobstructed 
views.  Hence,  an  alternative  is  to  use  motion  sensors  (for  example,  accelerometer,  rate 
gyroscope)  directly  attached  in  the  vicinity  of  center  of  mass  of  the  moving  platform  to 
acquire  acceleration,  velocity  and/or  displacement  as  a function  of  time.  For  this  case, 
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the  accuracy  of  the  location  of  the  sensor  could  be  the  question.  Hence  the  question  in 
this  section  is,  how  sensitive  is  the  calculated  Fa  to  errors  in  center  of  mass  locations 
estimates  where  the  identification  is  run  for  estimating  only  mass  and  inertia  parameters. 
Or,  with  limited  sensing,  use  other  measured  data  to  calculate  the  center  of  mass  motion. 
This  research  will  examine  the  sensor  location  effect  using  the  classical  sensitivity 
analysis. 

Using  equation  (2.55),  one  can  identify  the  dynamic  variables.  In  these  dynamic 
terms,  the  location  of  the  center  of  mass  is  chosen  for  checking  sensitivity.  Since  the 
location  of  the  center  of  mass  and  the  orientation  of  the  inertia  tensors  case  both  affect  the 
accuracy  of  the  identification  of  the  dynamics,  the  importance  of  finding  the  exact 
location  would  be  heightened.  False  measurement  of  the  center  of  mass  location  may 
cause  the  vagueness  of  the  results.  In  the  simulation,  the  starting  point  of  the  center  of 
the  mass  location  (G)  is  same  point  as  the  geometric  center  location  (C).  Point  Pj  is 
location  of  joint  i.  However  for  the  physical  system,  G will  actually  have  an  offset  ( ) 

from  C.  Now,  whether  the  offset  vector  r , plays  a role  for  the  sensitivity  of  the  force 

variation  is  in  question.  To  answer  this,  the  classical  dimensionless  sensitivity  equation 
(2.89)  is  used 

~ (A(Fa ).  )■ ) _ % change  in  (F  ).  (2 

r"  (A(u/)V(uJ,)  % change  in(rcg\ 

where  Faj  is  the  ith  component  of  the  actuator  force.  This  dimensionless  equation  was 
calculated  using  its  numerical  approximation.  The  s'r"'  represents  the  sensitivity  when 


the  offset  ( r ) exists. 

V eg  t 
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2.4.2.2  Trajectory  Selection  and  Conditions 

For  the  trajectory  condition,  one  needs  to  specify  which  direction  to  move.  For  this 
research,  sinusoidal  motions  are  selected  (Table  2.3).  Case  numbers  are  defined  for  the 
different  conditions  and  trajectory  directions.  Case  numbers  have  two  different 
trajectories.  For  pure  translation  case,  the  center  of  mass  location  moves  into  Xw 
direction  with  low  frequency  (0.5  Hz) 

2. 4.2.3  Sensitivity  Results 

The  geometric  center  location  of  special  6-6  PKM  is  [0, 0,-0. 35]  (meter).  From  this 
location,  a one  percent  change  of  the  maximum  distance  traveling  to  the  Xw  direction  is 

chosen  for  the  sensitivity  evaluation.  The  results  of  the  simulation  are  shown  in  Table 
2.4.  The  case  of  Rcg ^ 0 and  Beta_offset=0,  the  notation  x#y#z#  stands  for  the  directional 
effect.  The  ‘0’  means  no  change  in  that  direction  and  the  ‘ 1’  means  changes  were  made 
for  that  direction.  The  xOylzl  means  that  there  was  no  change  in  x direction  and  1 
percent  change  in  both  y and  z directions.  Note,  the  first  row  in  each  case  corresponds  to 
no  design  parameter  variation,  hence  is  always  zero. 

To  interpret  the  results  of  equation  (2.89)  for  each  case  shown  in  Table  2.4,  a value 
of  5.5  in  the  table  corresponds  to  a 5.5%  to  1%  ratio  of  change.  Any  value  greater  than 
one  is  considered  significant  since  the  design  parameter  variation  effect  is  amplified. 
When  Rcg  is  the  selected  design  parameter,  the  offset  effects  of  the  center  of  mass  are 
seen  to  be  greatest  in  the  terms  of  variations  in  the  design  parameters  X and  Z.  Even 
during  rotation,  directional  effects  were  of  the  same  trend  patterns. 

To  summarize  this  section’s  findings,  an  initial  study  of  the  sensitivity  of  the  sensor 
location  was  made.  The  location  offset  was  given  as  design  sensitivity  parameters.  The 
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findings  of  this  study  indicate  that  the  accuracy  of  the  center  of  mass  sensors  potentially 
has  a significant  effect  on  the  calculated  actuator  force  (Fa)  values  depending  on  how 

large  of  error  is  in  the  location  of  the  center  of  mass.  Hence,  simplification  of  the  system 
identification  algorithm  may  be  limited  in  terms  of  accuracy  of  results  desired.  Chapter  5 
presents  the  experimentally  obtained  identification  results  in  which  the  center  of  masses 
were  assumed  known.  Since  the  inertia  estimates  and  center  of  mass  locations  for  UF- 
PKM  have  been  identified  in  prior  research  and  the  results  readily  available  the 
assumptions  is  considered  not  too  conservative  nor  hinders  demonstrations  of  the  UKF 
identification  algorithms. 
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Table  2 1 Grouping  of  Terms  by  contributions  to  JFa 
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L XT  XT  J 

™,g 

k/  d 

All  remaining  me,  mb,  Ie,  and  Ib 
terms  noting.  TjJ,  [fJ  Tj  f, 

\ft  tt  f 

[f  T T T] 

L fft  grs  J 

m,g,  mh g 

k;  d 

Ffric  ,Tfrjc 

k:  tj] 

m a 

p c 

k.r  d 

F m/b 

[t  7 T t 

L cor  cor  J 

mZ£(m£xs3),  T(/£>o£ 

^4(®bxs3)’ 

, Zfs3  x v?)xs3 

k;  d 

Ftool 

[fr  T T\ 

mc 

Is 

e 

to 

c 

mc 

IBs 

mh 

m 

e 

o),  x [to,  x (T,s3  )JJ 

^kx^x’rj 
X I X I B^B 

x/£to£,m  (co/co £ 
to;jx[cDax(4s3)]] 

3 X [tot  X [tot  x J 

s,  X CO  x co  x“> 
S,  x CO  X co  x'V 

3 c c pc 

y*  3 

Ed 

,Exs3 

If  t t rf 

L umd  umd  J 

Humd  rTumd 

53 


Table  2.2  Motion  Trajectory  Cases  for  Section  2.2 


Motion  Trajectory  Cases 

Case  Number 

Trajectory 

Amplitude 

Radial 
Offset  Xw- 
Yw 

Direction 

1 

Sinusoidal  motion  in  X direction 

0.0254  m 

0m 

2 

Sinusoidal  motion  in  X direction 

0.0254  m 

0.0127  m 

3 

Sinusoidal  motion  in  Z direction 

0.0254  m 

0m 

4 

Sinusoidal  motion  in  Z direction 

0.0254  m 

0.0127  m 

5 

Sinusoidal  (0,  v|/)  rotational  motion 

5 degrees 

0m 

6 

Sinusoidal  (0,  v|/)  rotational  motion 

5 degrees 

0.0127  m 

[_ " 1 

Platform 

Motion 

Amplitudes 

Low  -> 
Medium  -> 
High  -> 

Velocity:  (0.0254)*7r  m/sec 
Velocity:  (0.0254)*27T  m/sec 
Velocity:  (0.0254)*47i  m/sec 

Acceleration:  (0.0254)*7i 
2 m/sec2 
Acceleration: 
(0.0254)*4ti2  m/sec2 
Acceleration: 

(0.0254)*  1 6tt2  m/sec2 

Frequency 

Low  -> 
Medium  -> 
High  -> 

n rad/sec 
2n  rad/sec 

4n  rad/sec 

Table  2.3  Motion  Trajectory  Cases  for  Section  2.4 


Case  Number 

Design  Parameter 
Variation 

Trajectory 

Amplitude 

1 

SF“ 

Pure  Translation,  X- 
direction 

0.0254  m 

2 

Pure  rotation 

5 degree 

Table  2.4  Table  of  Results  for  Section  2.4 
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Figure  2.1  Schematic  of  a Stewart-Gough  PKM 
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Figure  2.2  Physical  Diagram  of  Strut  i 
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Figure  2.3  Lumped  Mass  Parameter  Diagram  for  Strut  i 
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(a)  UF-SAMM,  Special  6-6  PKM  (side:SIDE=0.75) 


side  = 0.6096  m 
(b)  SNL  Hexel,  6-3  PKM  (side:SIDE=1.29) 


SIDE  = 0.4722  m ,0-2°14  m 


0.0382  m, 


side  = 0.5878  m 


0.1235  m SIDE  = 0.5276  m 

(c)  Ford  Research  Laboratory  PKM,  General  6-6  PKM  (side:SIDE=1.12) 


Figure  2 4 PKM  Configuration  Variants 
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Figure  2.6  Naming  Conventions 


Figure  2.7  Strut  Coordinate  System 
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Figure  2.8  Flowchart  for  the  UKF 
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Figure  2.9  Block  Diagram  of  Unscented  Transformation  [45] 
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Figure  2.10  Results  Defined  by  Equations  2.76  through  2.81  for  Trajectories  of  Table  2.2  (Low  Velocity  Cases) 

UF-PKM  with  side:SIDE=0.75  and  Ratio  mp/(6*me)=0.45 
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Figure  2.1 1 Results  Defined  by  Equations  2.76  through  2.81  for  Trajectories  of  Table  2.2  (High  Velocity  Cases) 

UF-PKM  with  side:SIDE=0.75  and  Ratio  mp/(6*me)=0.45 
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Figure  2.12  Results  Defined  by  Equations  2.76  through  2.81  for  Trajectories  of  Table  2.2  (High  Velocity  Cases) 

UF-PKM  with  side:SIDE=0.75  and  Ratio  mp/(6*me)=3.34 
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Figure  2.13  Results  Defined  by  Equations  2.76  through  2.81  for  Trajectories  of  Table  2.2  (High  Velocity  Cases) 

6-3-PKM  with  side:SIDE=1.29  and  Ratio  mp/(6*me)=0.45 
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Figure  2.14  Results  Defined  by  Equations  2.82  through  2.88  for  Trajectories  of  Table  2.2  (Low  Velocity  Cases) 

UF-PKM  with  side:SIDE=0.75  and  Ratio  mp/(6*me)=0.45 
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Figure  2.15  Results  Defined  by  Equations  2.82  through  2.88  for  Trajectories  of  Table  2.2  (High  Velocity  Cases) 

UF-PKM  with  side:SIDE=0.75  and  Ratio  mp/(6*me)=0.45 


CHAPTER  3 

IMPLEMEMTATION  OF  PARAMETER  IDENTIFICATION  METHODS 
This  chapter  illustrates  the  implementation  of  parameter  identification  methods  and 
provides  an  evaluation  of  their  performance.  The  inverted  double  pendulum  example 
was  selected  as  the  case  study  due  to  its  nonlinear  dynamics  and  the  availability  of  results 
by  others  in  the  literature.  Simulation  results  delineate  the  effectiveness  of  the  commonly 
used  least  squares  (LS)  method  and  the  recently  developed  unscented  Kalman  filter 
(UKF)  estimation  method  for  estimating  the  dynamic  parameters  of  nonlinear  systems. 
The  relative  advantages  and  disadvantages  of  one  method  versus  another  are  evaluated 
for  determining  each  method's  applicability  in  parameter  identification  of  parallel 
mechanisms.  This  comparison  is  augmented  with  the  results  from  the  literature  on  the 
use  of  the  extended  Kalman  filter  (EKF).  This  chapter  closes  with  the  decision  to  use 
UKF  for  parameter  identification  of  PKM  dynamics.  Implementation  results  for  the 
PKM  are  presented  in  Chapters  4 and  5. 

3.1  Pendulum  Setup  and  Equations  of  Motion 

An  inverted  double  pendulum  example  has  previously  been  used  for  comparing 
extended  Kalman  filter  (EKF)  and  unscented  Kalman  filter  (UKF)  theory  [38],  In  this 
dissertation,  the  model  and  conditions  of  references  [38,  45]  are  recreated  for 
demonstrating  the  implementation  of  the  UKF  as  well  as  adding  the  comparison  of  results 
obtained  with  the  LS  method.  To  make  the  simulation  study  more  realistic,  noise  was 

added  to  the  algorithm  s simulated  measured  states  (Xt);  thus,  emulating  measurement 


68 


69 


noise.  In  the  inverted  double  pendulum  example,  noise  was  not  added  to  the  measured 

output  ( dk ) since  the  simulated  system  had  zero  driving  forces  and  torques. 

The  inverted  double  pendulum  (Figure  3.1)  has  the  following  states  (x)  and 
parameters  (w) 


x=[x,  x,0x,6x,Q2,62\ 

(3.1) 

w =[ll,l2,m],m2,M] 

(3.2) 

: x : position  of  the  cart 

x : velocity  of  the  cart 

<9, : angular  position  of  the  bottom  pendulum,  link  1 (absolute) 

<9, : rate  of  change  in  angular  position  of  the  bottom  pendulum,  link  1 (absolute) 
d2 : angular  position  of  the  top  pendulum,  link  2 (absolute) 

02 : rate  of  change  in  angular  position  of  the  top  pendulum,  link  2 (absolute) 
h , m\ : length  and  mass  of  the  bottom  pendulum  (link  1),  respectively 
l2,m2:  length  and  mass  of  the  top  pendulum  (link  2),  respectively 
M : mass  of  the  cart 

The  mass  of  each  link  i is  assumed  to  be  uniformly  distributed,  where  the  center  of 

mass  is  at  A / . The  gravitational  force  is  assumed  to  be  in  the  vertical  direction  as  shown 

2 ' 

in  Figure  3.1.  The  dynamic  equations  of  motion  of  this  example  are  as  follows. 

(M  + ml  + m2 )x  + (—  m]  + m2)lxOx  cos <9,  + — m2l2 02  cos 02 
2 2 

= F + (-ml  + m2)ltO sin  6{  +-m2l202~  sin 

2 2 


(3.3) 
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(~ m\  +m2)lxxcos6x  + (~  + m2)lx 2 <9,  + -m2lxl202  cos(02  - 6 ,) 

1 1-2 
= (—  mx  + m2 )glx  sin 6X  +~m2lj202  sin(#2  -#,) 


(3.4) 


1 


1 


— m2xl2  cos<92  +-m2lj20]  cos(<92  -0])  + -m1l2  02 

1 1-2 
= —m2gl2  sin^ m2lxl2Qx  sin(6f  -6*,) 


(3.5) 


The  variable  (F)  is  the  control  input  force.  For  both  the  LS  and  UKF 
implementation,  a ‘freefalling’  test  setup  was  applied  to  obtain  response  data  for  the 
system  state  variables.  The  ‘freefalling’  test  denotes  generation  of  response  (system 
state)  data  without  any  controller  input  applied  to  the  system,  i.e.  only  external  load  is 
gravity.  In  a physical  implementation  of  an  identification  method,  this  data  would  be 
measured  experimentally.  Flowever,  to  emulate  the  identification  procedure,  the  system 
state  data  is  generated  using  the  ode45  solver  in  Matlab  for  the  system  as  modeled  in 
equations  (3.3)  thru  (3.5),  with  the  ‘true’  parameter  values  as  known.  The  data  collection 
time  interval  for  the  system  is  selected  to  be  0 second  to  20  seconds  where  F is  set  to  zero 
(no  control  input)  for  the  entire  simulation. 


3.2  The  Least  Squares  Method 

To  use  the  LS  method,  the  reference  residual  (e.  = {yk  - dk ). , eqn  (2.58))  must  be 
defined.  In  this  example,  the  dk  is  the  measured  control  input  force  and  y,  is  the 
corresponding  calculated  control  input  force  for  each  kth  iteration  using  the  assumed 
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theoretical  model  of  the  system  given  by  equations  (3.3)  thru  (3.5).  I.e., 


T*  = 


~F~ 

0 

= 

0 

k 

(M  + m,  + m2)x  + (—m,  +n 72 )/, <9,  cos#,  + — m2l202  cos02 
2 2 

-(— w,  + m2)l,0,2  sin  (9,  - — m2l202 2 sin#2 
2 2 

(- m, +m2)l,x  cos  0,  +(^  + m2V*fl  + k2/,/2#2  cos(#2  - #,) 


1 


1 • 2 

-(-"*,+  m2)gl,  sin  9,  --m2l,l202  sin(#2  - 0,) 

2 m2^2  c°s  ^ w2/,/2<9,  cos(#2  - 0, ) + ^ m2l20 


1 1-2 

- - m2gl2  sin  02  + - m2/,/26>2  sin(#2  - 0, ) 


(3.6) 


The  dk  is  set  to  [0  0 0]  (no  control  input  force)  for  all  time  in  the  ‘freefalling’  test.  The 
corresponding  JLSk  matrix  from  equation  (2.61)  is 


J IS.k  — 


Jr 


Jn  J» 
J 22  J 22 


J\A  J\S 

J 24  J 25 

J 24  J 22 


de, 

3ex 

de, 

de, 

de, 

dl, 

dl2 

dm 

dm2 

dM 

de2 

de2 

de2 

de2 

de2 

dl , 

dl2 

dm 

dm7 

dM 

Wko 

de2 

^2 

de. 

de. 

de2 

k 

dm, 

dm2 

dM 

(3.7) 


where 


1 


Ju  =(—ml  + m2  )0,  cos  0,  -(—m,  + m1  )0,  2 sin  0, 
2 2 


1 1-2 

J ,2  = —m20 2 cos 6*2  — m202  sin#2 
2 2 


1 ■■  l-, 

JX3  = x + — lx0]cosd]  ~—l\0\  sin#. 


(3.7a) 

(3.7b) 

(3.7c) 


-/  * 
2' 


Jl4  = x + 1,0,  cos0 , +—l202  cos02  -1,0,  sin0 , — 1202  sin#2 


(3.7d) 
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J \ 5 ~ * 


(3.7e) 


J 2i  = (— w,  + w2)xcos#,  +2(—  + m2 )/,6>!  h — m2l202  cos(<92  -0,) 


(3.7f) 


- (-^  w,  + m2 )g sin 9X  - — w2/2 /92 2 sin(<92  - <9, ) 

J2i  = ^m2lx02cos(02  -^1)-^-OT2/1^22sin(^2  -0,)  (3.7g) 

1 1 2 1 

J23  = -/,x cos#,  + -/,#,  --g/.sintf,  (3.7h) 


J24  = /,xcos<9,  + /,20,  + ~lj202  cos(02  -<9i)-^-/,/2(922  sin(02  -0,)-g/,  sin#. 


(3.7i) 

to 

II 

0 

(3.7j) 

1 1 . 2 

J 31  = —m2l20x  cos (#2  -#,)  + — m2l20x~  sin(#2  -#,) 

2 2 

(3.7k) 

1 1 2 

<7 32  = —m2x cos#2  H — m2/,6>1  cos(#2  -0x)  + — m2l202 
2 2 3 

(3.71) 

-^-w2gsin^2  +^-w2/,(9|2sin(6'2  -#,) 

^33=0  (3.7m) 

J34  = -/2xcos#2  + ~lj20x  cos(02 -0X)  + UX292  -X-gl2  sin#2  + ^-/,/2#,2  sin(#2  -#,) 

(3-7n) 

^35=0  (3.7o) 

As  stated  in  Section  2.3.1  (equation  2.64),  parameter  updates  are  calculated  each 
iteration  using 

Avr^  = (J LS,kJ LS,k  ) J LS,kek 


(3.8) 
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w*+i=w*+Awt  (3.9) 

where  e | is  the  kth  reference  residual  evaluated  at  w , — w . I.e., 

K wko  k ko  5 

eklt  = (yk~dXk„  (3-io) 


The  estimation  is  considered  to  have  converged  to  true  values  when  each  element  in 
ek  L,„  has  a value  less  than  the  acceptable  tolerance.  In  this  Chapter’s  simulations,  a 
tolerance  of  lxl O'7  was  used. 

3.3  The  Unscented  Kalman  Filter  Setting 

Since  the  UKF  uses  the  original  nonlinear  equations,  it  does  not  need  a 
linearization  procedure.  I.e.,  the  original  nonlinear  equations  of  motion,  (3.3)-(3.5),  with 
zero  control  input  yield 


(3.11) 


and 


(M  + mx  + m2  )x  + (—  m]  + m2 )/, <9,  cos  <9,  + — m2l2  02  cos  02 
7 7 


~(—m]  +m2)/16?,2  sin 6^,  - — m2l202  sin#. 


(t-w,  + m2 )/jXcos + w2)/,2(9|  + -m2lj2d2  cos(<92 -6>) 
^ 3 2 

-(|w,  +m2)glt  sin<9,  -^m2lj2022  sin(02 -0,) 


G(xk , w)  = 


— m2xl2  cost?,  +~m2lll20]  cos (<92  -0i)  + -m2l2  02 


1 . 1 

— m2gl2  sin  <9,  + — m2lj202  sin(02  -0,) 
2 2 


(3.12) 
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As  done  in  reference  [45],  all  the  parameters  were  initialized  to  1 .0.  The  other 
UKF  searching  values  were  selected  as  follows:  a = 0.0001 , p = 2,  k = 3-n,  where  for 

the  inverted  double  pendulum  example  the  dimension  of  w is  n=  5.  The  updating  of  R[ 
or  R[  uses  Robbins-Monro  Method  (see  section  2.3.2.3,  equation  2.75). 


3.4  Results 

As  can  be  seen  in  Table  3.1  for  the  LS  estimate  using  the  ‘without  noise’  model, 
the  /,  and  /2  parameter  values  converged  well  to  the  true  values,  however  the  remaining 
parameter  values  converged  only  to  within  1 5%  error  from  the  true  value.  A similar 
convergence  trend  can  be  found  for  the  extended  Kalman  filter  method  [45]  as  shown  in 
the  bottom  section  ot  Table  3.1.  For  the  LS  method,  a total  of  1225  iteration  steps  were 
required  to  reach  convergence  within  the  allowable  error  tolerance  (set  to  lxl O'7).  The 
total  CPU  time  of  this  run  was  4.23  seconds  with  a Pentium  4(1.8  GHz),  5 12  M RAM 
computer.  Figures  3.2  to  3.6  show  a plot  of  the  estimation  results  for  each  parameter  at 
each  iteration  step  of  the  LS  method  (dotted  lines).  The  dashed  line  denotes  the 
corresponding  true  parameter  value  and  the  solid  line  is  the  UKF  results. 

For  the  UKF  estimate  in  Table  3.1,  all  parameters  converged  to  true  value  within  a 
1 % error  range  in  2128  iterations  and  68.88  CPU  seconds  for  a tolerance  setting  of 
1x10  7 . However,  if  one  examines  Figures  3.2  to  3.6,  it  is  found  that  the  UKF 
converged  quickly  to  within  a close  range  of  the  error  tolerance  but  the  program  ran  for 
approximately  another  2000  iterations  and  took  68.88  CPU  seconds  before  reaching  the 
given  error  tolerance  ( 1 x 10  7 ).  Thus,  with  a reduction  in  tolerance  on  the  UKF 
convergence  criterion,  the  CPU  would  be  approximately  the  same  as  that  for  the  LS 
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method  (maintaining  original  tolerance).  As  seen  in  Table  3.1,  the  general  trend  in 
accuracy  of  parameter  estimation  in  terms  of  tolerance  for  both  LS  and  UKF  methods 
showed  a decrease  in  accuracy  as  the  tolerance  was  loosened  (e.g.  tolerance  lxl O'7  -» 

lxl  O'  resulted  in  accuracy  decreasing  by  factors  up  to  22  times  for  UKF,  yet  still  under 
1%  to  3.5%  error). 

Another  comparison  was  made  by  adding  noise  in  the  system.  To  emulate  the  real 
system,  random  noise  was  added  to  the  value  generated  from  the  data  generating  code 
described  in  Section  j.  1 thus  emulating  measurement  noise.  The  maximum  value  of  this 
random  noise  was  selected  to  be  less  than  1 % of  the  generated  value.  This 
correspondences  to  motion  sensors  with  error  specifications  as  follows:  for  x this 
corresponded  to  possible  maximum  error  of  3 centimeters  for  a 3 meters  range  of  motion; 
maximum  error  in  X is  0.02  meters  per  second  for  range  of  2.19  meters  per  second; 
maximum  q | is  0.36  radians  for  range  of  36.5  radians;  maximum  q is  0.18  radians  for 

range  of  18.2  radians  per  second;  maximum  X is  0.5  meters  per  second  squared  for  range 
of  5 1 . 1 meters  per  second  squared.  The  results  of  the  LS  and  UKF  methods  subject  to 
this  random  error  are  shown  on  Table  3.2  and  Figures  3.7  through  3.11.  The  UKF 
estimation  once  again  outperformed  the  LS  estimation  in  terms  of  accuracy  where  both 
had  the  same  number  of  iterations.  While  both  estimation  methods  have  merit,  the  LS 
method,  however  took  less  CPU  time  at  the  cost  of  requiring  a priori  linearization  of  the 
model  base.  For  some  parameters,  it  is  observed  that  the  measurement  noise  prevented 
premature  termination  of  the  LS  estimation  process,  yielding  slightly  improved  estimates. 

I he  EKF  outperforms  the  LS  but  it  too  is  limited  by  the  need  to  linearize  equations  of 
motion  before  implementation.  To  conclude,  although  the  UKF  has  a greater  overhead  in 
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CPU  time,  its  greater  convergence  accuracy  and  ease  of  implementation  for  nonlinear 


systems  (e.g.,  no  linearization  requirements)  makes  it  the  best  candidate  for  dynamic 


parameter  estimation  of  PKM’s. 


Table  3.1  Comparison  Between  UKF  and  LS  Without  Noise 
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Figure  3.1  Double  Inverted  Pendulum 


LS  vs  UKF  without  Noise 


Figure  3.2  Value  of  /, , Case  without  Noise  (LS  vs.  UKF) 
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Figure  3.3  Value  of  l2 , Case  without  Noise  (LS  vs.  UKF) 


LS  vs  UKF  without  Noise 


Figure  3.4  Value  of  m] , Case  without  Noise  (LS  vs.  UKF) 
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LS  vs  UKF  without  Noise 


Figure  3.5  Vaiue  of  m2 , Case  without  Noise  (LS  vs.  UKF) 
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Figure  3.6  Value  of  M , Case  without  Noise  (LS  vs.  UKF) 
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Figure  3.7  Value  of  /, , Case  with  Noise  (LS  vs.  UKF) 


Figure  3.8  Value  of  l2 , Case  with  Noise  (LS  vs.  UKF) 
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Figure  3.9  Value  of  mx , Case  with  Noise  (LS  vs.  UKF) 


Figure  3.10  Value  of  m2,  Case  with  Noise  (LS  vs.  UKF) 
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Figure  3.11  Value  of  M , Case  with  Noise  (LS  vs.  UKF) 


CHAPTER  4 

IMPLEMENTATION  OF  UKF  FOR  PARAMETER  IDENTIFICATION 

OF  PKM  DYNAMICS 

This  chapter  presents  a simulated  implementation  of  the  unscented  Kalman  filter 
(UKF)  for  parameter  identification  of  PKM  system  parameters.  The  UF  Special  6-6  PKM 
system  was  used  to  illustrate  the  performance  of  the  identification  algorithm  in  both 
simulated  (Chapter  4)  and  experimental  (Chapter  5)  implementations.  Implementation  of 
the  UKF  parameter  identification  algorithm  was  done  for  two  different  forms  of  the  PKM 
model  base:  force-based  model  (equation  2.55)  and  energy-based  model  (equation  2.56). 

The  key  difference  between  these  two  implementations  is  in  the  measured  data  i e d 

’ k 

desired  output,  and  Xk  measured  states.  The  force-based  model’s  measured  data  includes 
the  actuator  forces,  and  the  platform’s  linear  and  angular  positions,  velocities  and 
accelerations.  The  energy-based  model’s  measured  data  includes  the  actuator  forces, 
change  in  strut  lengths,  and  the  platform’s  linear  and  angular  position  and  velocities.  For 
the  experimental  implementation,  the  states  were  calculated  from  the  measured  strut 
lengths,  time,  and  numerical  derivatives  determined  from  filtered  link  length  measured 
data.  The  velocity  and  acceleration  data  can  be  either  measured  directly  or  determined 
from  differentiated  measured  position  data.  The  differentiation  process  can  further 
exacerbate  the  problem  of  noise  and/or  measurement  error  by  the  introduction  of  an 
additional  source  of  numerical  errors.  Since  the  data  acquired  from  experiments  always 
contains  some  form  of  noise  and/or  measurement  errors,  a comparative  analysis  of  the 
UKF  parameter  identification  algorithm’s  performance  for  the  two  different  model-bases 
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was  conducted  to  determine  if  there  is  any  benefit  in  using  the  energy-based  model 
approach  other  than  reduced  modeling  effort  in  deriving  the  theoretical  model.  The 
anticipated  benefit  is  reduced  noise  effects  due  to  lower  order  derivatives  and/or  the 
number  of  measurements  required  in  the  energy-based  method.  If  the  results  show  a 
significant  enhancement  in  performance,  then  a greater  potential  existed  for  the  use  of  the 
UKF  identification  method  for  online  estimation.  For  an  industrial  setting,  online 
estimation  can  be  the  critical  factor  that  potentially  enables  the  engineer  to  further 
improve  a machine  s performance  by  reducing  the  operational  time  and  machining  errors. 

4.1  PKM  Simulations  (Force-Based) 

Referring  to  the  PKM  model  presented  in  Chapter  2,  the  measured  data  in  the 
torce-based  model  are  the  actuator  forces,  and  the  platform's  linear  and  angular  positions, 
velocities  and  accelerations.  To  obtain  the  emulated  measured  actuator  forces  for  a given 
state  response,  a simulation  data-generating  program  of  the  model  (equation  2.55)  was 
used  in  which  the  emulated  measured  states  were  entered  into  the  program  as  given.  The 
emulated  measured  states  are  (refer  to  equations  (2.4)  through  (2.8)) 

x*  =lqk  >q„  J (4.1) 

The  emulated  measured  (desired)  output  data  is 

— (Fa  )mamircJ  (4.2) 

and  the  calculated  output  is 

y*=G(xk,wk)  (4.3) 

Again,  it  is  noted  that  tor  implementation  purposes,  the  center  of  mass  locations  were 
assumed  known,  i.e.,  not  included  in  the  set  of  parameters  to  be  identified.  In  this 
dissertation,  the  CAD-based  values  that  were  determined  by  reference  [48]  are  used.  In 
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the  analysis  and  implementation,  it  was  also  assumed  that  all  struts  are  of  modular  design 
and  have  the  same  inertia  parameters,  thus  reducing  the  number  of  parameters  in  the 


identification  set.  Hence,  the  parameters  identified  ( Wk ) in  the  implementations  are 


w 


= [me^b,mp,IxxbJyyh,IzzbJxxpJ 


,/  ,/  I 

5 yye  5 zzc  J* 


(4.4) 


As  stated  above,  the  data-generating  program  consists  of  the  system  model  defined 


in  Section  2.1  where  the  'true'  values  for  the  inertia  parameters  (CAD-based  values  in 
this  dissertation)  were  used  for  generating  the  emulated  measured  data.  In  this  chapter, 
two  different  tests  were  conducted:  freefalling  test  and  actively  controlled  strut  test.  The 
freefalling  test  was  used  to  generate  emulated  measured  data  as  was  done  for  the  example 
problem  in  Chapter  3.  The  actively  controlled  strut  test  emulates  an  experimental 
implementation  in  which  the  strut  forces  are  actively  controlled.  The  freefalling  test 
might  be  considered  more  useful  in  cases  where  the  strut  force  data  is  very  noisy  or 
contains  significant  errors.  After  acquiring  the  emulated  measured  data  from  the  data- 
generating  program  for  each  of  these  tests,  the  data  was  fed  into  the  UKF  algorithm  for 
estimating  the  mass  parameters.  To  further  build  confidence  in  the  UKF  parameter 
estimation  method  for  PKMs,  friction  is  introduced  into  the  system  and  the  results  are 
compared  to  the  frictionless  case.  As  a representative  model  of  friction  found  in 
mechanical  systems,  the  modified  friction  model  [46]  was  selected  and  is  defined  as 
follows. 


F 


frict 


Fried,  1)  • L + Fric(i,2) 
• (2  • Fried, 2) / 0.004) • Z. 
Fric(i,\)  ■ L - Fric(i, 2) 


for  L > 0.002 

for  - 0.002  <L  <0.002  <4-5) 

for  L < -0.002 


where  Fried, \)  = viscous  friction  coefficient  for  ith  strut 
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Fric(i,2)  = Coulomb  friction  coefficient  for  ith  strut 


L,  - velocity  of  individual  strut  (Unit:  m/s) 


Table  4.1  provides  the  estimated  coefficients  for  the  friction  model  obtain  using  single 
strut  tests  of  reference  [37],  Figure  4.1  shows  the  vertical  strut  friction  testbed  and  a 
representative  plot  of  the  results  obtained  for  strut  1 [37], 


The  differences  between  the  actual  implementation  of  force-based  model  and 
energy-based  model  in  UKF  parameter  identification  are  the  model-base  and  the 
measured  data.  All  other  assumptions  hold  for  both  energy-based  and  force-based 
implementations.  In  the  implementation  cases  presented,  it  is  noted  that  the  energy-based 
model  results  were  tor  only  the  actively  controlled  strut  tests.  The  energy-based  model’s 
data-generating  program  is  the  same  as  the  one  used  for  force-based  model  with  an 
additional  emulated  measured  data  set  of  (A L.  )k  for  each  strut  i and  time  step  k.  Hence, 
the  emulated  measured  states  for  equation  (2.56)  are 

x*  =(qkT»q/)r  (4.6) 

and  the  emulated  measured  actuator  ((desired)  measured  output)  force  data  is 


where  the  parameters  identified  (Wk  ) are  the  same  set  as  used  for  force-based. 


4.2  PKM  Simulations  (Energy-Based  Model) 


(4.8) 


(4.7) 


xxb  5 yyb  5 


(4.9) 
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4.3  Simulation  Setup 

The  same  physical  setup  and  conditions  for  both  the  force-based  modeling  and 
energy-based  model  implementations  were  used  for  each  simulation  case.  For  the 
treefalling  test,  the  movable  platform  starts  at  the  'home'  position  (initial  condition).  The 
'home  position  is  defined  as  the  PKM  configuration  for  which  the  mass  center  of  the 
movable  platform  is  located  at  [x  = 0.00,  y = 0.00,  z - - 0.35]  meters  from  the  inertial 
frame  { W}  as  defined  in  Section  2. 1 .3. 1 . The  physical  dimensions  of  the  PKM  are  as 
defined  in  Figure  2.4  (a). 

For  the  actively  controlled  strut  test,  the  movable  platform  starts  at  the  ‘home’ 
position  (initial  condition).  The  mass  center  of  the  movable  platform  is  then  commanded 
to  move  sinusoidally  about  the  'home'  position  in  the  x-direction  of  the  inertial  frame 
with  amplitude  of  one  inch  (case  1).  It  should  be  noted  that  the  maximum  width  of  UF- 
PKM  workspace  is  approximately  4 inches  in  the  Xw  and  Yw  plane.  For  case  2,  the 
movable  platform  moves  sinusoidally  with  the  same  amplitude  about  the  ‘home’  position 
in  the  z-direction  of  the  inertial  frame.  For  case  3,  the  center  of  mass  of  the  movable 
platform  moves  sinusoidally  about  the  ‘home’  position  in  the  positive/negative  x and 
positive/negative  y direction.  The  frequency  of  the  movement  is  4*7c  Hz  (approximately 
45  % of  the  maximum  speed  capability  of  the  UF  PKM).  The  initial  conditions  of  all 
motions  are  contained  in  Table  4.2.  The  trajectories  (cases)  used  in  the  simulations  and 
their  parameter  values  are  specified  in  Table  4.3 

Again,  note  it  has  been  assumed  that  all  struts  have  the  same  inertia  parameters  and 
that  the  center  of  mass  locations  and  kinematic  parameters  are  known.  The  rule  of 
notation  used  in  the  simulations  are  as  follows: 
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• All  headers  and  titles  in  the  tables  and  figures  relating  to  force-based  model  use  the 
notation  type:  FTYPE1TYPE2  TYPE3 

o For  FTYPE1 

■ FA  denotes  ‘active  controlled  strut  tests’ 

■ FF  for  ‘freefalling  tests’ 

■ Note  that  for  FF  type,  there  is  no  TYPE2. 

o For  TYPE2 

■ X for  x-directional  movement,  Z for  z-directional  movement,  and 
XY  for  xy-directional  movement 

o For  TYPE3 

■ W indicates  ‘with  friction’  and  WO  represents  ‘without  friction’ 

• All  notations  relating  to  energy-based  model  follow  the  same  rule  as  force-based 
model  with  the  exception  of  using  E instead  of  F for  denoting  the  TYPE1  case. 

4.4  Simulation  Results  (Force-Based,  Without  Friction) 

Table  4.4  and  4.5  show  the  simulation  results  for  the  force-based  model,  without 
friction.  For  mass  parameters,  compared  to  the  true  value,  the  x-direction  movement 
(FA_X_WO)  and  the  xy  direction  movement  (FA  XY  WO)  yield  the  best  estimates. 

The  estimated  mass  values  converged  within  0.04%  of  its  true  value.  The  estimated  mass 
parameters  for  z direction  movement  (FA_Z_WO)  and  freefalling  case  (FF  WO) 
however  converged  to  within  16.5%  and  1 1.4%,  respectively.  Table  4.5  shows  the  value 
of  convergence  for  inertia  parameters.  In  all  cases,  the  movable  platform  parameter 
estimates  (Ixxp,  Iyyp,  Izzp)  had  better  convergence  than  strut  parameter  estimates  (Ixxe,  Iyye, 
Izze,  fixb-  Iyyb,  Fzzb)-  The  maximum  percentage  difference  for  platform  parameters  is  0.8%, 
however,  the  maximum  percentage  difference  for  strut  parameters  is  26.2%.  The  largest 
percentage  difference  can  be  seen  on  Iyye  This  could  be  observed  in  two  ways.  One,  the 
numerical  value  assigned  in  the  simulation  was  very  small  compared  to  other  parameter 
values  (i.e.,  the  largest  parameter  in  the  simulation  was  7.7214,  which  is  roughly  142,000 
times  greater  than  Iyye  value  (inertia  about  longitudinal  axis  of  ball  spline)).  And/or,  the 
selected  test  trajectories  (data)  do  not  yield  sufficient  dynamic  excitation  for 
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identification  ot  the  Iyye  parameter.  The  dynamic  excitation  effects  were  explored  further 
by  randomizing  or  spreading  the  response  data  inputs  ( x, , d, ).  These  effects  will  be 
discussed  further  in  Section  4.9. 

Figures  4.2-4.15  are  plots  of  the  results  of  the  force-based,  with  and  without 
friction.  These  plots  are  grouped  into  four  categories:  me  (Figures  4.2  and  4.3),  mb 
(Figures  4.4  and  4.5),  mp  (Figures  4.6  and  4.7),  and  all  mass  parameters  (Figures  4.8- 
4.1 1).  The  plots  in  Figures  4.3,  4.5,  4.7  and  4.8-4. 1 1 (without  friction  case)  show  a 
parameter  convergence  tendency  towards  their  true  values  for  the  actively  controlled  strut 
cases  (FA_X,  FA  XY).  For  me  case,  the  results  of  both  x directional  movement  and  xy 
directional  movement  converged  faster  and  with  greater  accuracy  than  the  other  two  (z 
and  freefalling)  trajectories.  This  tendency  can  also  be  found  in  other  mass  cases  (mb, 
mp).  For  freefalling  test,  it  shows  generally  poorer  convergence  accuracy  for  all  cases. 
This  can  be  attributed  to  the  limitation  of  the  strut  and  the  workspace  of  the  UF  Special  6- 
6 PKM  on  its  range  of  motion  as  well  as  the  PKM  configuration  approaching  a geometric 
condition  in  which  struts  are  essentially  parallel.  That  is,  since  the  range  of  strut  lengths 
is  limited,  the  freefalling’  movement  has  only  a limited  amount  of  data  since  the  motion 
occurred  too  fast  to  capture  in  detail.  This  data  was  not  sufficient  to  run  the  unscented 
Kalman  filter  algorithm  for  extended  period  without  either  reusing  data  sets  or  truncating 
the  iterations.  In  this  dissertation  the  iterations  were  truncated. 

4.5  Simulation  Results  (Force-Based,  With  Friction) 

1 able  4.6  and  4.7  show  the  values  of  convergence  for  the  force-based  model,  with 
friction  results.  For  the  mass  properties,  the  estimated  values  are  in  error  from  the  true 
values  by  less  than  3%  (x-direction),  2.6%  (z-direction)  and  16.2%  (xy-direction).  The 
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convergence  error  between  the  estimated  value  and  true  value  in  xy-direction  exhibited  a 
large  increase  compared  to  the  without  friction  case.  The  results  of  the  simulation  for  the 
‘with  friction'  case  are  shown  for  me  case  (Figure  4.2),  mb  case  (Figure  4.4),  mp  case 
(Figure  4.6),  and  all  mass  parameters  (Figures  4.12  through  4.15).  For  me  case,  z- 
direction  case  is  the  fastest  and  closest  to  the  true  value.  The  x-direction  cases  and  xy- 
direction  converge  close  enough  to  be  considered  good,  but  show  slower  convergence 
than  the  z-direction  trajectory.  Similar  trends  were  found  on  mb  and  mp  case. 

4.6  Simulation  Results  (Energy-Based,  Without  Friction) 

C omparing  only  the  cases  without  friction,  the  energy-based  model  results  are 
shown  in  Tables  4.8  and  4.9  for  the  estimated  mass  and  inertia  parameters.  The  plots 
showing  parameter  convergence  tendencies  are  in  Figures  4.16  thru  4.2 1 . The  resulting 
percentage  differences  between  the  true  and  estimated  values  all  remained  under  7.3  % 
for  the  mass  properties.  In  general,  the  convergence  performance  rate  and  accuracy  using 
energy-based  model  case  was  worse  than  that  achieved  using  force-based  model  case 
with  exception  of  the  mb  estimate  for  the  FA_Z  _WO  case.  It  is  also  observed  in 
comparing  energy-based  results  without  friction  with  the  force-based  results  with  friction, 
the  convergence  accuracies  exhibited  to  a large  degree  the  same  order  of  magnitude  in 
percentage  error.  One  notable  difference  (seen  from  results  in  Tables  4.4  and  4.8)  is  that 
the  energy-based  model  achieved  much  better  convergence  accuracy,  yielding  a lower  the 
percentage  difference  between  true  and  estimated  parameter  value  of  Iyye.  The  percentage 
differences  between  true  and  estimated  for  the  other  inertia  parameters  are  of  either 
similar  or  higher  magnitudes  then  that  obtained  using  the  force-based  model.  The 
energy-based  model  case  was  not  done  using  the  freefalling  and  the  with  friction  tests. 
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hence  assuming  similar  degrade  in  convergence  accuracy  would  result  as  that  found  for 
the  force-based  model  results. 

The  convergence  tendencies  are  shown  on  Figures  4.16  (me),  4.17  (mb),  4.18  (mp). 
For  me  case,  x and  xy  directional  movement  shows  the  similar  trend  of  convergence. 

This  trend  can  also  be  seen  in  mb  and  mp  cases.  The  z-directional  movement  shows  more 
of  a linear  tendency  in  convergence. 

4.7  Simulation  Results  (Force-Based,  Without  Friction,  Original/Random/Spread) 

The  k measured  data  ( , d k ) consists  of  the  system  states  and  actuator  forces 
measured  at  a given  instant  in  time.  In  theory,  given  a perfect  model-base  and  noise-free 
data,  one  could  use  a single  set  of  measured  data  (x^,d4 ) (i.e.,  (x<  ,dt)=(xA  ,dt  , ) for 

all  k)  and  run  the  UKF  algorithm  until  wt  converges  to  the  true  estimate.  Also,  since 
Wk  's  an  estimate  of  non-time  varying  parameters  of  the  system,  then  any  kth  set  of 
measured  data  (xt  ,dt ) should  yield  the  same  w,  estimate.  However,  depending  on  the 

values  ot  the  measured  states,  it  is  possible  that  some  Wk  parameters  would  not  be 

identifiable  because  at  the  given  instant  in  time  their  contribution  to  the  dynamics 
response  is  zero  or  essentially  negligible.  Thus,  in  practice,  at  each  kth  iteration  a new  set 
ot  measured  data  (xk , d* ) is  used  so  as  to  ensure  the  effects  of  all  identifiable  parameters 
are  captured  in  the  estimation  process.  The  use  of  multiple  data  sets  (x;. , dA ) is  also 

necessary  when  these  sets  are  obtained  from  experiments,  that  in  general  contain  noise 
and  contributions  from  the  unmodeled  dynamics. 

There  are  infinity  possible  data  sets  ( xt , dt. ) from  which  one  can  select.  However 

there  is  no  assurance,  that  the  selected  data  sets  will  yield  a — » ‘true’  values,  the 
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question  arises  as  to  which  collection  of  data  sets  (xt  ,dA ) should  one  choose.  That  is, 
what  motion  trajectory  should  the  PKM  be  driven  to  follow  so  as  to  yield  a sufficiently 
dynamically  exciting  response  from  which  to  extract  (xk , dt ).  In  search  of  a guarantee, 
researchers  have  developed  various  methodologies  for  generating  sufficiently  exciting 
trajectories  [25],  As  in  reference  [45],  van  der  Merwe  et  al.  used  a free-fall  approach. 
Per  the  study  conducted  in  section  2.4,  it  was  determined  that  sinusoidal  trajectories 
about  the  “home  ’ position  provide  sufficiently  dynamically  exciting  data.  However,  to 
further  substantiate  the  capabilities  and  limitations  in  using  UKF  and  to  explore  the 
dynamic  excitation  effects,  the  force-based  model  UKF  method  was  implemented  where 
the  response  data  time  incremental  sets  were  either  fed  into  the  algorithm  in  a random 
time  order  or  simply  in  an  increased  spread  in  the  time  increment  of  the  selected  data. 
Figure  4.22  illustrates  the  difference  in  the  ordering  of  the  data  sets  for  one  variable 
input.  The  results  in  the  prior  sections  used  response  time  ordered  ( x; . , A k ) data  sets  as 

depicted  by  the  original  case.  To  restate,  since  the  time  increments  of  the  previous 
simulations  from  one  time  incremental  step  to  the  next  were  very  small  (0.001  second), 
the  dynamic  variation  of  the  response  data  was  limited  as  well  throughout  the  estimation 
process.  Thus,  randomizing  or  spreading  the  (xt  ,dt)  data  sequence  was  used  as  a 

means  tor  validating  that  the  prior  results  had  sufficient  dynamic  excitation  and  to  gauge 
how  robust  the  UKF  method.  The  term  ‘random’  means  that  the  sequence  of  the 
( xk  5 d* ) data  sets  given  to  the  estimation  algorithm  is  randomly  chosen  and  the  term 
‘spread’  represents  that  the  sequence  of  ( xt , dt ) data  sets  given  to  the  estimation 
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algorithm  is  chosen,  e.g.,  from  every  10  incremental  steps  (i.e.  instead  of  using  data 
number  1,  2,  3,  ...  , using  1,  1 1, 21,  31,  ...). 

For  the  force-based  model  without  friction,  Tables  4.10  and  4.1 1 show  the  results 
using  randomized  data.  The  percentage  difference  in  x and  xy  direction  shows  good 
convergence.  However,  similar  to  the  case  using  sequential  data,  the  z-directional 
movement  and  freefalling  tests  do  not  yield  high  convergence.  Tables  4.12  and  4.13 
contain  the  results  for  ‘spread'  case.  These  results  show  better  performance  in  the 
percentage  difference  compared  to  the  ‘random’  case.  Thus,  the  x and  xy  directional 
movements  show  superiority  in  convergence  performance  where  the  z direction  exhibited 
the  worst  performance. 

The  plots  that  display  the  convergence  of  the  ‘random’  and  ‘spread’  case  are  in 
Figures  4.23  thru  4.42  where  plots  of  the  original  estimation  results  (from  Figures  4.8 
through  4.1 1)  are  at  times  indistinguishable  from  the  spread  plots.  That  is,  all  figures 
relating  to  me,  mb,  and  mp,  respectively,  in  x direction  and  xy  direction  case  present 
similar  results.  Data  from  the  original  (sequential)  set  and  the  spread  case  are 
comparable  to  each  other,  and  the  random  case  exhibits  a slower  convergence  tendency. 
The  z-direction  case  and  treefalling  cases  do  not  have  good  convergence  when  using  the 
data  in  the  random  nor  the  spread  case,  again  experiencing  the  same  limitations  as  cited 
earlier. 

4.8  Comparison  Between  Force-  and  Energy-Based  with  Noisy  Data 

As  an  additional  comparison  between  the  force-based  model  and  the  energy-based 
model  cases  without  friction,  noise  was  added  into  the  UKF  estimation.  The  1 % noise 
(per  definition  in  Chapter  3)  was  added  to  the  UKF  measured  data  ( xk ) of  the  states.  The 
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force-based  model  case  was  also  ran  with  5 % noise.  As  a representative  case,  the  results 
for  the  xy-directional  trajectory  without  friction  simulation  case  is  chosen,  shown  in 
Tables  4. 14  through  4. 17.  These  tables  show  that  as  the  noise  increases  the  difference 
between  the  estimates  and  the  CAD-based  ‘true’  values  increases  as  expected.  Figures 
4.43  and  4.44  also  support  this  observation.  When  comparing  Tables  4. 14  and  4.16,  the 
mass  parameters  are  estimated  within  5.4  % for  force-based  model  and  37.8  % for 
energy-based  model  when  the  1%  noise  is  added.  This  result  shows  that  the  force-based 
model  outperforms  the  energy-based  approach  and  hence  is  the  preferred  model-base  to 
be  used  for  the  UKF.  Additional  convergence  results  for  the  other  mass  parameters  are 
provided  in  Appendix  C. 

4.9  Summary  of  Results 

Summarizing  this  Chapter  provided  comparisons  on  the  selection  of  the  model-base 
for  UKF  parameter  identification  in  terms  of  estimation  convergence  accuracy  and  rate 
trends,  for  cases  without  friction,  with  different  motion  trajectories,  and  with  and  without 
noise.  It  was  found  that  the  x-directional  case  exhibited  the  lowest  percentage  difference 
between  estimated  and  true  parameter  values.  The  xy-direction  cases  had  similar  trends 
for  the  x-direction  cases;  however,  when  friction  was  added  to  the  system,  the  results 
were  not  as  reliable  as  those  found  for  the  x-direction  case.  The  convergence  tendency 
for  the  z-directional  movement  and  the  ‘freefalling’  test  were  both  poor,  in  general. 
Mostly,  ‘freefalling’  tests  suffered  from  the  limited  range  of  motion  in  which  to  capture  a 
sufficient  amount  of  data  as  well  as  the  directional  symmetry  of  the  strut  axes  (also  seen 
for  the  z-directional  movement). 

For  the  z-direction  cases,  energy-based  model  is  preferred  to  the  force-based  model 
due  to  a lower  percentage  difference  between  the  estimated  and  true  values  of  the 
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relatively  small  parameters.  However,  for  x-direction  and  xy-direction  case,  force-based 
model  is  superior  to  the  energy-based  model  in  terms  of  accuracy  of  convergence. 

Hence,  a significant  improvement  using  energy-based  method  was  not  generally  realized. 
For  data  set  types  (original/random/spread),  the  results  indicate  that  original  data  or 
spread  data  usage  is  the  preferred  for  obtaining  better  estimation  results,  as  well  as  it 
requires  less  data  processing  effort. 
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Table  4.4  Mass  Value  (Units  - kilogram)-  Force-Based,  Without  Friction 
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Table  4.6  Mass  Value  (Units  - kilogram)-  Force-Based,  With  Friction 
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Table  4.8  Mass  Value  (Units  - kilogram)-  Energy-Based,  Without  Friction 
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Table  4.10  Mass  Value  (Units  kilogram)—  Force-Based,  Without  Friction  Randomized 
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Table  4.12  Mass  Value  (Units  - kilogram)-  Force-Based,  Without  Friction  Spread 
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Table  4.14  Mass  Values  with  Noisy  Data  (Units  - kilogram)-  Force-Based,  Without  Friction  (0,  1 and  5 % Noise) 
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(a)  The  Vertical  Strut  Friction  Test  Stand 


(b)  Strut  1 Friction  Data 


Figure  4.1  Single  Strut  Friction  Test  [37] 


lTle  Value  (estimated  value)  - Force  with  Friction 


Figure  4.2  me  Value  (Force  with  Friction) 
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me  Value  (estimated  value)  - Force  without  Friction 


Figure  4.3  me  Value  (Force  without  Friction) 


mbValue  (estimated  value)  - Force  with  Friction 


Figure  4.4  mb  Value  (Force  with  Friction) 
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m^Value  (estimated  value)  - Force  without  Friction 
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Figure  4.5  mb  Value  (Force  without  Friction) 


m pValue  (estimated  value)  - Force  with  Friction 


Figure  4.6  mp  Value  (Force  with  Friction) 
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nipValue  (estimated  value)  - Force  without  Friction 
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Figure  4.7  mp  Value  (Force  without  Friction) 
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Figure  4.8  Mass  Values  (X-direction  without  Friction) 
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Mass  Parameters  (estimated  value)  - FA  XY  WO 
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Figure  4.9  Mass  Values  (XY-direction  without  Friction) 
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Figure  4.10  Mass  Values  (Z-direction  without  Friction) 
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Figure  4.1 1 Mass  Values  (FF  case  without  Friction) 
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Mass  Parameters  (estimated  value)  - FA  X W 


, 1 1 1 1 1 

mb 

Ime 

imb 

Pr 

— True  Value 

i 

1 1 

3 

o 

_J 1 I 

/ 

1 1 1 l I i i 

Figure  4.12  Mass  Values  (X-direction  with  Friction) 
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Figure  4.13  Mass  Values  (XY-direction  with  Friction) 
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Figure  4.14  Mass  Values  (Z-direction  with  Friction) 
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Mass  Parameters  (estimated  value)  - FF  W 
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Figure  4.15  Mass  Values  (FF  case  with  Friction) 


Hie  Value  (estimated  value)  - Energy  without  Friction 


Figure  4.16  me  Value  (Energy  without  Friction) 
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mi,  Value  (estimated  value)  - Energy  without  Friction 


Figure  4.17  mb  Value  (Energy  without  Friction) 


mp  Value  (estimated  value)  - Energy  without  Friction 


Figure  4.18  mp  Value  (Energy  without  Friction) 
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Figure  4.19  Mass  Values  (X-direction:  Energy-Based  without  Friction) 


Figure  4.20  Mass  Values  (XY-direction:  Energy-Based  without  Friction) 
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Mass  Parameters  (estimated  value)  - case  EA  X WO 
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Figure  4.21  Mass  Values  (Z-direction:  Energy-Based  without  Friction) 


(5)  (1)  (7)  (3)  (2)  (5)  (4) 

Response  time  increment 
(interaction  # k) 

Figure  4.22  Data  Set  Collection  Type  (Original,  Spreading  and  Randomizing) 
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me  Value  (estimated  value)  - FA  X case 


Figure  4.23  me  Value  (FA  X WO  - Original/Random/Spread) 


ITlt)  Value  (estimated  value)  - FA  X case 


Figure  4.24  mb  Value  (FA  X WO  - Original/Random/Spread) 
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nip  Value  (estimated  value)  - FA  X case 


Figure  4.25  mp  Value  (FA  X WO  - Original/Random/Spread) 


me  Value  (estimated  value)  - FA  Z case 


Figure  4.26  me  Value  (FA  Z WO  — Original/Random/Spread) 
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mi, Value  (estimated  value)  - FA  Z case 


Figure  4.27  mb  Value  (FA  Z WO  - Original/Random/Spread) 


lTlpValue  (estimated  value)  - FA  Z case 


Figure  4.28  mp  Value  (FA  Z WO  - Original/Random/Spread) 
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me  Value  (estimated  value)  - FA  XY  case 


Figure  4.29  me  Value  (FA  XY  WO  - Original/Random/Spread) 


ni(,  Value  (estimated  value)  - FA  XY  case 


Figure  4.30  mb  Value  (FA  XY  WO  - Original/Random/Spread) 
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IllpValue  (estimated  value)  - FA  XY  case 


Figure  4.31  mp  Value  (FA  XY  WO  - Original/Random/Spread) 


me  Value  (estimated  value)  - FF  case 


Figure  4.32  me  Value  (FF  WO  - Original/Random/Spread) 
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1T1|,  Value  (estimated  value)  - FF  case 


Figure  4.33  mb  Value  (FF  WO  - Original/Random/Spread) 


trip  Value  (estimated  value)  - FF  case 


Figure  4.34  mp  Value  (FF  WO  - Original/Random/Spread) 
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Mass  Parameters  (estimated  value)  - FA  X WO 


y 

i 1 1 , — — 

— n\ 
mb 

8 

mb 

nip 

True  Value 

7 

i ' 

- 

- 

6 

A 

- 

E 

2 

o)  5 
o 

i \ 
i ^ 

i ^ 

mp 

5 

4 

- 

3 

- 

2 

me 

1 

1 j — 1 i i 

-1 1 1 1 i 

10  20  30  40  50  60 


Iteration 


Figure  4.35  Mass  Value  (X-direction,  Random) 
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Figure  4.36  Mass  Value  (XY-direction,  Random) 
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Mass  Parameters  (estimated  value)  - FA  Z WO 
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Figure  4.38  Mass  Value  (FF  case.  Random) 
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Figure  4.39  Mass  Value  (X-direction,  Spread) 
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Figure  4.40  Mass  Value  (XY-direction,  Spread) 
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Figure  4.41  Mass  Value  (Z-direction,  Spread) 
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Figure  4.42  Mass  Value  (FF  case,  Spread) 
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me  Value  (estimated  value)  - Force  without  Friction  - with  Noise  (0,  1,5  percent) 
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Figure  4.43  me  Value  (0,  1 and  5 % Noise,  Force-Based  Model) 
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Figure  4.44  me  Value  (0  and  1 % Noise,  Energy-Based  Model) 


CHAPTER  5 

PKM  EXPERIMENT  RESULTS 


This  chapter  presents  the  experimental  setups  and  the  estimated  inertial  parameter 
results  for  the  UF-PKM.  Since  the  energy-base  method  did  not  perform  as  well  as 
expected,  only  the  force-based  UKF  approach  is  used  in  this  dissertation  for  the 
experimental  validation  study.  An  ultimate  goal  for  this  dissertation  is  to  determine  how 
well  the  estimation  program  performs  in  a real  PKM  application  with  machine  and 
instrumentation  limitations.  Since  the  strut  lengths  and  individual  strut  forces  are  the 
most  readily  available  data  that  can  be  acquired  from  the  machine  controller  and  sensors, 
the  accuracy  of  this  data  is  a critical  for  extracting  other  system  states  such  as  strut 
velocities  and  accelerations.  Shamblin  [37]  noted  that  the  encoder  reading  for  acquiring 
strut  length  data  and  current  reading  for  acquiring  the  strut  force  contained  significant 
noise  resulting  in  measurement  inaccuracies.  This  noise  effect  was  further  amplified  as 
numerical  approximations  were  used  in  determining  the  strut  velocities  and  accelerations. 
Therefore,  to  implement  the  UKF  identification  procedure  presented  in  Chapters  2 
through  4,  sources  of  error  and  limitations  are  expected  such  as  those  listed  below. 

• Signal  Noise 

• Kinematic  Errors 

• Measurement  Errors 

• Controller  Parameter  and  Modeling  Errors 

• Directional  Control  Problem 

• Following  Error 

• Poor  Encoder/  Current  Readings 

• Limited  Sensing  and  Actuation 
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Hence,  in  the  experimental  implementation  of  UKF,  corresponding  uncertainties  and 
reconciliation  techniques  had  to  be  taken  into  consideration  throughout  the  analysis  and 
identification  of  the  system  parameters.  As  discussed  in  reference  [47]  and  chapter  2,  the 
effects  of  limited  sensing  were  considered  in  regards  to  assuming  the  geometric  center  of 
the  platform  could  be  assumed  the  same  as  the  center  of  mass.  Theoretical  results 
indicate  that  errors  in  this  approximation  would  result  in  errors  in  the  resulting  parameter 
estimation.  However,  due  to  the  simplicity  and  the  full  knowledge  of  the  PKM’s 
components  properties,  the  implementation  of  the  UKF  parameter  identification 
algorithm  was  done  assuming  all  centers  of  masses  per  reference  [48].  Assuming  ideal 
kinematics,  the  limited  sensing  issues  (e.g.,  non-direct  measurement  of  platform  motion) 
were  circumvented  using  kinematics  and  numerical  derivatives  to  determine  the  platform 
position,  velocity,  and  acceleration. 

5.1  Experimental  Setup 

The  UF  PKM  (refer  to  Figures  2.4  (a)  and  5.1)  is  actuated  by  Aerotech  BM200  3- 
phase  brushless  motors  with  encoder  feedback.  The  motors  are  powered  by  Aerotech 
BA20  amplifiers  using  pulse-width  modulation  to  provide  current  control.  The  motion 
control  software  is  the  open-architecture,  PC-based  Enhanced  Machine  Controller, 
developed  by  the  National  Institute  of  Standard  and  Technology  (NIST),  and  runs  under 
Real-Time  Linux  on  a Pentium  II  PC.  The  interface  between  the  amplifiers,  motors,  and 
software  is  a Servo-to-Go  motion  controller  card.  At  the  present  level  of  controller 
development,  only  translational  motions  of  the  platform  are  achievable. 

In  the  experiments,  the  platform  was  moved  in  coordinated  motion  by  executing  a 
single  G-code  command  or  a series  of  commands  from  a script  file.  The  built-in  logging 
feature  was  used  to  record  strut  lengths  and  motor  current  values.  From  the  data,  actual 
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platform  position  for  each  time  step  was  obtained  using  an  iterative  forward  kinematics 
solver  [36].  The  platform  velocity  and  acceleration  were  calculated  using  the  PKM 
kinematics  and  the  strut  length,  velocities  and  acceleration.  The  strut  velocities  and 
accelerations  were  determined  numerically  applying  the  4-point  central  difference 
smoothing  equations  ([37],  equations  (5.1)  and  (5.2))  to  filtered  position  data. 


(2.1 1))  and  At  is  the  incremental  time  step. 

The  motor  current  readings  were  converted  to  the  measured  Fa  using  the  motor 
torque  constant,  and  transmission  ratio  of  the  ballscrew.  The  friction  forces  along  the 
strut  axes  were  estimated  using  strut  velocity  data  with  a Coulomb  and  viscous  friction 
model.  In  prior  research  on  the  UF-PKM  [37],  the  model’s  friction  coefficients  were 
obtained  experimentally  through  system  identification  of  each  individual  strut. 


Numerous  sources  of  error  presented  difficulties  in  experimentally  validating  the 
identification  of  the  dynamic  parameters.  Much  of  the  error  was  indirectly  caused  by  the 
limitations  of  the  servo  controller.  The  PKM’s  controller  uses  a PID  algorithm  to 
generate  the  control  effort  from  the  encoder  measured  errors.  In  tuning  the  system’s 
controller,  the  derivative  term  was  found  necessary  for  minimizing  oscillations  that 
would  occur  if  only  P (proportional)  or  PI  (proportional  - integral)  control  was  used.  The 
derivative  term  also  amplified  the  high  frequency  vibrations  that  arose  from  encoder 


(5.1) 


(5.2) 


where  (Lk  ).  is  the  absolute  strut  length  for  the  ith  strut  at  kth  iteration  step  (equation 


5.2  Challenges  to  Experimental  Validation  and  Reconciliation 
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noise.  The  result  is  that  the  actual  machine  motion  is  characterized  by  small  amplitude, 
high  frequency  positional  oscillations  in  addition  to  the  smooth  commanded  motion. 

In  addition  to  the  undesired  oscillations  are  the  ‘false’  movements  caused  by  finite 
encoder  resolution.  As  the  motor  rotates,  the  position  is  read  as  a series  of  step  moves  by 
the  encoder.  When  the  true  position  is  between  encoder  pulses,  the  position  as  read  by 
the  encoder  is  seen  to  oscillate  between  the  two  pulses,  introducing  oscillations  that  are 
not  actually  occurring.  Naturally,  differentiating  position  data  to  obtain  velocity  and 
especially,  acceleration,  amplifies  both  the  real  and  false  oscillations.  This  noise 
propagation  is  such  that  acceleration  caused  solely  by  encoder  noise  is  of  similar 
magnitude  to  the  acceleration  of  the  commanded  motion.  For  this  study,  an  8th  order 
Butterworth  filter  was  utilized  to  eliminate  positional  noise  of  frequencies  greater  than 
1 0Hz,  as  the  commanded  machine  movements  were  of  much  lower  frequency.  This 
alleviated  much  of  the  problem,  but  noise  at  less  than  10Hz  still  existed  and  became 
indistinguishable  from  actual  motion. 

Errors  due  to  noise  are  also  manifested  in  empirically  determined  constants  such  as 
friction  coefficients.  During  a friction  test,  although  the  commanded  velocity  is  constant, 
the  actual  velocity  and  the  corresponding  control  effort  are  oscillatory.  This  contributes 
to  random  error  when  measuring  the  electrical  current  required  to  produce  a constant  strut 
velocity.  Measurement  of  the  motor  torque  constant  (used  for  converting  current 
measurements  into  strut  forces)  was  similarly  effected. 

Another  source  of  error  is  unmodeled  transmission  dynamics.  For  analysis 
purposes,  the  transmission  of  motor  rotation  to  linear  strut  movement  is  assumed  constant 
and  rigid.  The  actual  machine  uses  a flexible  rubber  coupling  between  the  motor  shaft 
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and  ballscrew,  acting  as  a spring-damper  system.  Thus,  oscillations  occurring  in  the 
motor  would  not  be  completely  transmitted  to  the  platform  motion,  as  is  assumed  in  the 
simulation.  Additionally,  the  ballscrew  itself  acts  as  a rotational  spring,  though  its  effect 
on  error  is  much  less  than  the  flexible  coupling. 

Compounding  the  previous  errors  in  the  friction  and  torque  constants  is  the 
relatively  small  magnitude  of  measured  dynamic  forces  in  relation  to  friction  forces. 

Thus,  to  reduce  friction  within  the  system,  3-to-l  precision  gearboxes  were  removed  from 
the  original  UF-PKM  setup.  Removing  the  gearbox  also  had  the  effect  of  reducing  motor 
inertia  effects  (Fm/b)-  Furthermore,  because  of  the  dominance  of  the  friction  with  the  UF- 
PKM,  these  and  other  unmodeled  dynamics  effects  were  assumed  negligible  in  the 
implementation  of  the  UKF  identification  method.  Due  to  time  limitations,  only  two  sets 
of  single  strut  friction  coefficients  were  updated  with  new  experimental  results  (strut  1 
and  4).  The  single  strut  friction  tests  were  performed  using  a horizontal  single  strut 
configuration.  This  eliminated  the  effects  from  gravity.  The  resultant  estimated 
coefficients  for  the  friction  model  using  single  strut  tests  are  shown  on  Table  4. 1 . Using 
these  results  compared  to  reference  [37],  reference  [37]  ’s  estimates  were  scaled  to 
provide  new  estimates  for  the  remaining  struts,  thus  accounting  for  the  change  in  the  gear 
ratios  with  the  new  experimental  setup. 

5.3  Experiment  Results  (Without  Gearbox) 

As  stated  previously,  the  experimental  data  for  this  dissertation  is  acquired  from 
experiments  conducted  with  the  gearbox  removed  from  the  PKM  system.  In  the  fully 
assembled  UF-PKM  experiments,  coordinated  x,  xy,  and  z-directional  movements  were 
used  to  generate  the  response  and  output  data  sets.  The  PKM  platform’s  direction  and 
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speed  commanded  were  the  same  as  those  used  for  the  simulated  analysis  as  shown  Table 
4.3,  excluding  the  freefalling  test. 

5.3.1  Friction  Extraction  and  Estimation  Stages 

Since  the  friction  coefficients  (both  viscous  and  Coulomb)  also  need  to  be 
estimated  using  the  UKF  algorithm,  the  total  number  of  parameters  increased  to  24  (3 
masses,  9 principal  inertias  and  1 2 friction  coefficients)  The 

w*  = [ w*  ’ ^™C(U),  Fric{\,2),  Fric(2,\),  Fric(2,2),  Fric{\\),  Fric(3,2),  ^ 

Fric{  4,1),  Fric(  4,2),  Fric{  5,1),  Fric(  5,2),  Fric(  6,1),  Fric(  6,2)] 

where  it  is  assumed  that  all  struts  have  identical  mass  and  inertia  properties.  The 

principal  axes  of  inertia  for  each  strut  are  defined  along  the  {s2,s3,s23}.  axes,  and  the 

center  of  mass  locations  are  those  estimated  in  reference  [48].  The  parameter 
identification  was  first  executed  for  the  entire  set  of  24  parameters.  Using  the  following 
UKF  initial  estimates  to  start  the  estimation  process  (Reference  [37]  single  strut  friction 
estimates  plus  within  90%  of  CAD-based  mass  and  inertia),  the  parameter  estimation  did 
not  converge.  This  led  to  implementing  a three-stage  identification  approach. 

• Stage  1.  Conducting  an  estimation  of  friction  followed  by  inertia  parameter 
estimation.  Stage  1 has  two  steps: 

o Stage  1-A.  The  initial  guess  for  the  friction  estimation  is  made  (e.g.,  the 
prior  strut  estimates  [37]  and/or  the  author’s  single  strut  results).  At  this 
stage,  the  CAD-based  inertia  parameters  were  used  as  the  ‘true’  inertia 
values.  All  three  trajectories  (x,  z,  and  xy-direction)  were  tested  for 
estimating  the  friction  parameters. 
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W*L,,  = [^c(U),  Fried!),  Fric(2,\),  Fric(2,2), 

Fric(  3,1),  Fric{  3,2),  Fric(  4,1),  Fric(4,2),  (5-4) 

F>/c(5,l),  Fric(5,2),  Fric(  6,1),  Fn'c(6,2)f 


o Stage  1-B.  After  obtaining  the  estimated  friction  values  for  each  trajectory, 
the  inertia  parameters  then  are  estimated  with  the  estimated  friction 
parameters  from  Stage  1-A  used  as  the  ‘true’  friction  values.  All  three 
trajectories  are  conducted  as  well  to  estimate  the  inertia  parameters. 


wkL,-*  =K  ,mh,mp,Ix 


' I yyh  5 L 


• ’ I xxp  ’ I yyp  ’ 1 zip  ' I 


xxe  ’ I yye  ' ^ ~ 


■ ]/  (5-5) 


Note,  that  the  estimated  friction  parameters  from  FA  X (Stage  1-A)  were  used 
as  the  ‘true’  friction  values  in  estimating  the  inertia  parameters  in  FA  X. 
Also,  the  estimated  friction  parameters  from  FA_Z  (Stage  1 -A)  were 
used  as  the  ‘true’  friction  values  in  estimating  the  inertia  parameters  in 
FA  Z. 


• Stage  2.  Conducting  an  estimation  of  friction  followed  by  inertia  parameter 
estimation  with  using  the  data  given  in  Stage  1-A.  Stage  2 has  two  steps: 
o Stage  2-A.  The  initial  guess  for  the  friction  estimation  is  the  values 
from  Stage  1 -A.  At  this  stage,  the  CAD-based  inertia  parameters  were 
used  for  the  ‘true’  inertia  values.  All  three  trajectories  (x,  z,  and  xy- 
direction)  were  tested  for  estimating  the  friction  parameters. 

w*L2-,=w*L,-,  (5.6) 

o Stage  2-B.  After  obtaining  the  estimated  friction  values  for  each 

trajectory,  the  inertia  parameters  then  are  estimated  with  applying  the 
estimated  friction  parameters  from  Stage  2-A  used  as  the  ‘true’  friction 
values.  All  three  trajectories  are  conducted  as  well  to  estimate  the 
inertia  parameters. 
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w 


^ \ Stage  2 -B 


w 


^ \ Stage  \-B 


(5.7) 


Note,  that  the  estimated  friction  parameters  from  FA_X  (Stage  2-A)  were  used 
as  the  friction  values  to  estimate  the  inertia  parameters  in  FA_X.  Also 
the  estimated  friction  parameters  from  FA_Z  (Stage  2-A)  were  used  as 
the  friction  values  to  estimate  the  inertia  parameters  in  FA_Z. 

• Stage  3.  Both  the  friction  and  inertia  parameters  were  estimated 
simultaneously,  i.e., 

w*  = [ w * lorapto-4  ’ Fn'c(U),  Fric{\,2),  Fric(2,l),  Fric(2,2), 

Fric(3,l),  Fric(3,2),  Fric(4,\),  Fric(4,2),  (5-8) 

Fric(  5,1),  Fric(  5,2),  Fric(  6,1),  Fric(  6,2)] 

The  data  used  in  this  stage  as  the  initial  estimate  is  the  combined  values  from 

Stage  2-A  and  Stage  2-B.  For  example,  the  estimated  friction  values  from 

Stage  2-A  and  the  estimated  inertia  values  from  Stage  2-B  in  FAX  were 

combined  into  one  set  and  used  as  the  initial  parameter  values  to  estimate  both 

the  friction  and  the  inertia  parameters  in  FA_X,  FA  Z,  and  FA  XY  runs  for 

Stage  3. 

5.3.2  Experimental  Results  with  Friction  Estimates  (Stages  1 and  2) 

Table  5.1  shows  the  estimated  friction  value  using  the  friction  estimation  program 
(Stage  1-A).  This  estimation  uses  as  its  initial  value  those  shown  in  Table  4.1  and  the 
CAD  values  of  reference  [48]  as  the  true  values.  As  can  be  seen,  in  x-directional 
movement  (denoted  by  FA_X),  the  estimated  Coulomb  friction  value  ( Fric{i, 2) , refer  to 
equation  (4.1))  was  different  from  the  expected  single  strut  estimated  value  with  ~5  to 
250  % error.  For  the  estimated  viscous  friction  value  ( Fric(i,  1) , refer  to  equation  (4.1)), 
strut  1 and  6 coefficients  have  quite  different  values  from  the  single  strut  estimated  values 
with  ~54  and  32,  respectively,  % error  considering  all  trajectory  cases.  For  z-directional 
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movement  (denoted  by  FA_Z),  the  estimated  Coulomb  friction  values  are  different  in 
most  struts  with  ~5  to  210  % error.  The  viscous  frictions  for  strut  1,  2,  4,  and  6 have  an 
error  from  the  expected  value  with  range  of  ~36  to  85  %.  In  xy-directional  movement 
(denoted  by  FA_XY),  the  viscous  values  are  worse  than  the  previous  cases  (up  to  1 23  % 
error).  The  Coulomb  friction  differences  from  the  expected  values  are  in  the  range  of 
~15  to  300%.  Using  these  newly  estimated  values,  the  inertial  parameter  estimation  was 
performed  as  Stage  1-B  in  estimation.  The  results  for  Stage  1-B  are  shown  on  Table  5.2. 
The  estimated  mass  parameters  have  errors  ranging  from  ~2  to  19  % compared  to  the 
expected  CAD  based  estimates.  The  inertial  parameters  have  errors,  which  range  from 
0.7  to  13  % except  Iyye  value,  which  had  the  same  convergence  difficulty  as  in  simulation 
(Chapter  4). 

Since  the  percent  difference  in  the  friction  estimations  (Stage  1-A)  compared  to 
the  values  in  Table  4.1  is  too  large  and  the  difference  between  each  trajectory  values  are 
also  too  large,  the  next  stage  (Stage  2)  of  estimation  is  required  for  identifying  which 
friction  parameter  values  are  more  reliable.  The  results  are  shown  on  Tables  5.3  and  5.4. 
Note,  that  for  x-directional  movement;  the  previously  estimated  friction  parameters  of 
Stage  1 -A  from  the  x-directional  movement  case  were  used  as  the  expected  true  values 
for  the  friction  estimation,  and  similarly  for  z and  xy-directional  movements.  The 
estimated  friction  results  for  x-directional  movement  now  have  an  error  range  of  0 to  4 % 
from  previously  estimated  friction  values.  For  z-directional  movement,  the  error  range 
for  friction  was  1 5 %.  Generally,  an  error  less  than  2.6%  was  found  for  xy  motion  with 
exception  of  one  coefficient  having  ~23%  error.  With  using  the  resultant  friction  values 
from  Stage  2-A  (Table  5.3),  inertia  parameters  are  estimated  (Stage  2-B,  Table  5.4). 
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Again,  the  CAD-based  values  are  assumed  the  ‘true’  values  as  basis  for  comparison  of 
estimated  results.  In  addition,  the  Stage  1-B  values  were  compared  with  the  estimated 
values  to  determine  if  any  convergence  to  the  new  estimates  is  occurring.  For  mass  and 
inertia  parameter  values,  the  error  ranges  from  the  CAD  values  are  relatively  close  to  the 
previous  estimates  in  Table  5.2  (within  13%  to  18%).  Comparing  Stage  1-B  estimated 
values  with  Stage  2-B  estimates  indicates  that  there  is  possible  convergence  starting  to  be 
reached  for  each  trajectory  direction,  refer  to  (%  diff  Stage  1-B)  values  in  Table  5.4. 

5.3.3  Experimental  Results  with  Different  Initial  Values  (Stage  3) 

While  the  friction  estimates  are  converging  to  values  in  each  respective  direction, 
there  does  not  appear  to  be  a convergence  to  a single  set  of  parameter  values  for  friction 
regardless  of  motion  trajectory.  Hence,  Stage  3 is  a search  for  a single  set  of  friction 
parameter  values  using  full  parameter  estimation.  This  involves  using  each  of  the 
estimates  in  friction  of  Stage  2-A  and  the  estimated  inertia  values  from  Stage  2-B  for 
each  trajectory  as  the  initial  values  for  the  other  trajectories  in  Stage  3 (e.g.,  the  use 
values  from  FA  X as  the  initial  value  for  the  FA_XY  trajectory  estimation).  The 
rationale  for  this  is  one  to  avoid  estimates  to  have  reached  a local  minimum  as  well  as  to 
determine  an  estimate  of  a single  set  of  parameters.  If  the  program  falls  into  a local 
minimum  inside  the  UKF,  then  the  convergence  could  be  different  from  what  it  is 
supposed  to  be.  The  results  using  x-directional  value  are  shown  on  Tables  5.5  and  5.6. 
Use  of  other  directional  movement  (xy  and  z)  estimates  as  initial  values  were  found  to 
not  yield  convergence  in  the  estimation,  thus,  the  corresponding  initial  values  (estimates 
of  Stage  2-A  and  Stage  2-B)  are  not  estimates  of  the  ‘true’  values.  The  observation  of  the 
results  in  x-directional  movement  estimation  case  is  that  the  viscous  friction  values  all 
converged  within  a 0 to  22  % range,  but  the  Coulomb  friction  value  went  up  to  0 to  1 86 
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% for  the  FA_XY  trajectory.  The  inertial  parameters  converged  within  5 to  26  % from 
the  CAD  values.  However,  for  FA_X  and  FA_Z  trajectory  case,  when  comparing  to  the 
Stage  2-B  value,  the  inertia  parameters  converged  within  5 %.  For  the  FA_XY  trajectory 
case,  the  inertia  parameters  converged  within  39  %.  Figure  5.2  thru  5.16  show  the 
convergence  tendency  for  mass  parameters  and  friction  parameters.  Most  of  figures  have 
similar  tendencies.  The  xy-directional  movement  values  converged  to  much  higher 
values  than  the  z-directional  movement  values,  where  the  x-directional  movement  values 
converge  to  the  expected  values  (estimates  from  Stage  2 in  x-direction).  The  z- 
directional  movement  values  converge  much  closer  to  the  expected  values  than  the  xy- 
directional  movement  case.  For  Figure  5.12,  it  can  be  seen  that  although  the  intermediate 
estimates  diverged  from  expected  at  some  stage  of  iteration,  its  value  returned  quickly  to 
an  estimate,  converging  to  within  -8.67  % of  expected. 

5.4  Summary  of  Results 

For  ‘without  gearbox’  case,  several  observations  are  made  for  estimating  the 
parameters.  For  friction  case,  the  direct  measurement  from  the  individual  struts  using  the 
single  strut  tests  does  not  capture  the  spherical  and  Hooke  joint  friction  values.  This  can 
be  delineated  from  the  estimation  results  for  the  different  input  trajectories.  E.g.,  the  X- 
and  XY-  motions  move  the  struts  over  different  joint  ranges  while  the  Z-motion  yields 
joint  motion  closer  to  the  X-motion  case.  Selecting  as  initial  values  the  estimates  of  prior 
estimations  determined  along  different  motion  trajectories  quickly  narrowed  the  possible 
results  towards  finding  the  true  values  of  the  system.  The  inertia  parameters  exhibited  a 
much  faster  convergence  to  a ‘true’  value.  To  conclude,  the  single  strut  estimates  have 
not  fully  captured  the  strut  friction.  Thus,  the  new  friction  estimates  are  considered  more 
accurate.  The  inertia  estimates  showed  very  good,  accurate  convergence  especially  when 
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considering  the  fact  that  friction  forces  dominated  the  data.  However,  the  differences 
between  friction  estimates  per  the  different  trajectories  still  needs  further  resolve  as  to  the 
determination  of  the  true  rotational  friction  differences  of  each  spherical  and  Hooke 
joint’s  value.  Hence,  the  accurate  modeling  for  spherical  and  Hooke  joints  are  highly 
recommended  for  better  explanation  of  these  possible  phenomena  and  has  been  left  for 
future  work. 

To  conclude,  the  recommended  identification  methodology  for  actual 
implementation  on  PKM’s  is  the  multistage  iterative  approach  presented  in  section  5.3.1 
using  the  force-based  UKF  model  algorithm  of  Chapter  2.  This  is  particularly  true  for 
systems  with  a large  amount  of  friction.  The  following  highlights  some  of  the  critical 
issues. 

• First,  one  should  not  ignore  the  importance  of  the  initial  guess  for  the 
friction  parameters  (the  usual  dominance  found  in  most  ballscrew  driven 
machines).  Thus,  the  tests  with  different  trajectories  (i.e.  x,  z,  and  xy- 
directional  movement)  must  be  conducted  for  finding  the  possible  friction 
parameter  values.  (Stages  1 and  2) 

• Since  the  estimation  program  could  fall  into  the  local  minimum  for  the 
highly  nonlinear  systems,  the  estimated  friction  parameter  values  obtained 
from  the  previous  step  must  be  tested  again  with  other  trajectories,  (e.g.  the 
estimated  friction  parameters  from  x-directional  movement  test  must  be 
performed  with  xy-  and  z-directional  movement  test. 
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• The  final  friction  parameters,  that  are  found  to  show  estimateion 

convergence  in  all  other  trajectories,  should  be  taken  into  consideration  as 
the  closer  value  to  the  true  value. 


Table  5.1  Friction  Estimation  with  CAD  Based  Inertia  Assumed  True  Values  (Stage  1-A) 
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Table  5.3  Estimated  Friction  Value  Using  Previously  Estimated  Friction  Value  from  Stage  1-A  as  Initial  Value  fStape  ?-A3 
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Table  5.4  Estimated  Inertial  Value  Using  Previously  Estimated  Friction  Value  from  Stage  2-A  (Stage  2-B) 
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Table  5.5  Estimated  Friction  Value  Using  Previously  Estimated  Parameters  as  Initial  Value  (Stage  3) 
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Table  5.6  Estimated  Inertial  Value  Using  Previously  Estimated  Value  (x-directional  movement  from  Tables  5.3  and  5.4  (Stages  2-A 
and  2-B))  as  Initial  Value  (Stage  3)  
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Figure  5.1  Single  Strut  Testing 


mg  Value  (estimated  value) 


Figure  5.2  Value  of  me  (Experimental  Result) 
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mb  Value  (estimated  value) 


Figure  5.3  Value  of  mb  (Experimental  Result) 


Figure  5.4  Value  of  mp  (Experimental  Result) 


148 


Figure  5.5  Friction  (1,1)  Parameter  Value  (Experimental  Result) 


Figure  5.6  Friction  (1,2)  Parameter  Value  (Experimental  Result) 
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Figure  5.7  Friction  (2,1)  Parameter  Value  (Experimental  Result) 


Figure  5.8  Friction  (2,2)  Parameter  Value  (Experimental  Result) 
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Figure  5.9  Friction  (3,1)  Parameter  Value  (Experimental  Result) 


Figure  5.10  Friction  (3,2)  Parameter  Value  (Experimental  Result) 
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Figure  5.11  Friction  (4,1)  Parameter  Value  (Experimental  Result) 


Friction(4,2)  Parameters  (estimated  value) 
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Figure  5.12  Friction  (4,2)  Parameter  Value  (Experimental  Result) 
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Figure  5.13  Friction  (5,1)  Parameter  Value  (Experimental  Result) 


Figure  5.14  Friction  (5,2)  Parameter  Value  (Experimental  Result) 
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Figure  5.15  Friction  (6,1)  Parameter  Value  (Experimental  Result) 


Figure  5.16  Friction  (6,2)  Parameter  Value  (Experimental  Result) 


CHAPTER  6 
CONCLUSION 

To  date,  very  little  research  on  the  dynamics  of  parallel  kinematic  mechanisms 
(PKMs)  exists.  Of  the  work  available  in  the  literature,  most  have  been  for  generating  the 
equations  of  motion,  inverse  dynamics,  and  structural  dynamics.  Insight  to  the  structure 
and  properties  of  PKM  dynamics,  effects  of  strut  inertias,  dynamics  interactions  with 
various  kinematics  error  sources,  formulation  of  dynamic  performance  indices,  and 
dynamic  balancing  for  PKMs  have  only  recently  been  provided.  To  address  the 
increasing  awareness  of  the  dynamics  needs,  which  are  brought  to  the  forefront  by  the 
above  research  or  lack  thereof,  the  research  of  this  dissertation  established  a dynamic 
parameter  identification  approach  for  PKMs.  Thus,  this  enables  engineers  to  easily 
acquire  the  values  of  the  parameters  necessary  for  obtaining  a valid  dynamic  model  of  the 
system. 

To  choose  the  method  used  in  this  research,  several  identification  methods  were 
considered.  The  main  comparisons  were  made  between  least  squares  method,  extended 
Kalman  filter  and  unscented  Kalman  filter.  Using  the  inverted  double  pendulum  as  a 
sample  test  system  in  the  comparison  study,  it  was  found  that  the  UKF  performance  was 
overall  superior  to  the  others.  Thus,  UKF  was  selected  and  implemented  for  parameter 
identification  of  PKM  and  tested  via  simulation  study  and  validated  with  experimental 
implementation  on  the  UF  PKM. 

The  simulations  studies  for  the  UF  Special  6-6  PKM  had  two  different  categories: 
with  friction  and  without  friction.  Prior  work  indicated  a high  degree  of  friction  within 
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the  UF-PKM  system.  The  trajectories  were  chosen  to  exhibit  x,  z,  and  xy  direction 
relative  to  the  inertia  frame.  Other  comparisons  were  made  involving  different  tests  for 
generating  response  data,  e.g.,  actively  controlled  strut  test  and  ‘freefalling’  test. 

The  outcomes  of  the  simulations  were  that  the  estimated  mass  values  generally 
converged  well  within  1 percent  and  for  inertia  (excluding  Iyye)  values  converged  within  1 
to  12  percent  for  systems  ‘without-friction’  test.  When  the  values  of  the  inertia  were  very 
small  compared  to  other  parameter  values  (e.g.,  Iyye),  the  convergence  percentage  of  the 
inertia  was  only  within  26  percent.  For  ‘with-friction’  test,  the  convergence  trend  was 
similar  to  the  ‘without-friction’  test.  The  other  results  show  that  ‘freefalling’  test  did  not 
yield  improved  identification  over  the  use  of  the  ‘actively  controlled  strut’  test.  This  lack 
of  expected  improvement  is  in  part  attributed  to  the  limitation  of  UF  Special  6-6  PKM 
range  of  motion  as  well  as  the  geometric  effects  of  struts  approaching  a parallel 
configuration.  The  ‘freefalling’  test  needs  more  data  if  it  is  to  be  used  with  this 
methodology  is  to  be  used.  In  addition,  force-based  and  energy-based  model  methods 
were  compared  to  determine  if  there  is  any  benefit  other  than  reduced  model  effort  in 
deriving  the  theoretical  model.  A significant  improvement  using  energy-based  method 
was  not  generally  realized.  As  an  additional  comparison  between  the  force-based  model 
and  the  energy-based  model  cases  without  friction,  noise  was  added  into  the  UKF 
estimation.  The  mass  parameters  are  estimated  within  9.2  % for  force-based  model  and 
80.5  % for  energy-based  model  when  the  1%  noise  is  added.  This  result  shows  that  the 
force-based  model  is  the  highly  preferred  model-base  to  be  used  for  the  UKF. 

The  experimental  study  for  the  UF  Special  6-6  PKM  had  three  stages  to  be 
conducted.  In  Stage  1,  separate  parameter  sets  were  estimated  sequentially  (i.e.,  friction 
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parameter  estimation  followed  by  inertia  parameter  estimation  using  initial  friction 
parameter  CAD  based  inertia  estimates).  In  Stage  2,  again  involves  separate  parameter 
sets  sequentially  estimated  using  Stage  1 friction  parameter  estimates  as  initial  values  for 
the  UKF.  Stages  1 and  2 are  compared  for  determining  potential  convergence.  In  Stage 
3,  the  results  of  Stage  2 are  used  as  initial  values  in  UKF  for  full  parameter  estimation 
(both  friction  and  inertia  parameter).  The  initial  parameters  using  each  directional 
estimates  of  Stage  2 were  found  either  to  converge  or  not.  Thus,  the  outcome  allows 
quickly  narrows  the  possible  results  by  eliminating  the  non-convergence  case’s  initial 
values  as  possible  estimates  of  the  True’  value.  The  experimental  implementation’s 
inertial  parameter  estimation  was  considered  to  exhibit  high  degree  accuracy,  especially 
in  the  presence  of  high  friction  conditions  and  a large  degree  of  uncertainties  in  kinematic 
parameter. 

In  conclusion,  the  development  of  an  enhanced  dynamic  identification 
methodology  for  PKM  was  established  and  implemented.  The  simulated  results  show  the 
unscented  Kalman  filter  successfully  applied  for  identifying  the  system  parameters  of 
PKMs.  Since  the  unscented  Kalman  filter  uses  the  original  nonlinear  equations,  the 
procedure  for  obtaining  a linearized  model  is  not  necessary.  This  merit  corresponds  to 
the  significant  reduction  in  modeling  effort  required  to  generate  the  basic  algorithm  used 
for  the  system  identification.  In  an  industrial  setting,  since  a simple  method  to  be  used 
would  be  required,  the  unscented  Kalman  filter  would  be  the  best  choice.  The 
recommended  identification  methodology  for  actual  implementation  on  PKM  is  the 
multistage  iterative  approach  using  the  force-based  UKF  model.  This  is  particularly  true 
for  systems  with  a large  amount  friction.  The  highlights  of  some  of  the  critical  issues 
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were  presented  in  the  summary  of  Chapter  5.  In  brief  they  were  various  trajectories  must 
be  utilized  so  as  to  obtain  better  estimates  of  parameters  especially  in  the  presence  of 
friction  and  to  avoid  estimation  process  falling  into  local  minimum. 

The  overall  impact  to  technology  is  provision  of  a parameter  identification  method 
that  yields  improved  rigid  body  dynamics  models  of  parallel  kinematic  mechanisms. 

This  will  lead  to  the  ability  to  introduce  advanced  controllers  that  exploit  the  model 
dynamics  to  improve  system  performance.  In  addition,  with  some  refinement,  the 
presented  method  sets  the  foundation  for  formulating  an  online  tuning  algorithm. 

Other  future  variants  to  the  research  worthy  of  exploration  is: 

• Additional  testing  to  include  identification  of  the  center  of  mass  position  in 
the  UKF  experimental  implementation. 

• Dual  estimation  testing  (parameter  and  state  estimation). 

• Use  of  alternate  strut  force  sensing  (e.g.,  strain  gage  testing). 

• On-line  adaptive  tuning. 

• More  rigorous  approach  for  generating  dynamically  exciting  trajectories. 
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APPENDIX  A 

THE  Z-X-Z  EULER  ANGLE  DEFINITION  AND  [U]  MATRIX  FORMATION 
The  z-x-z  Euler  angles  are  defined  as  the  ordered  sequence  of  rotations  of 
coordinate  systems  with  rotations  about  z-axis,  rotated  x-axis,  and  rotated  z-axis.  The 
procedure  of  the  z-x-z  Euler  angle  rotations  are  as  follows  (as  shown  in  Figure  A.l). 

1 . The  base  coordinate  system  (e  coordinate  system,  same  as  {F}  coordinate  system 
in  Chapter  2)  is  rotated  to  f coordinate  system  with  angle  <j>  with  respect  to  z-axis 
( e3 ).  (refer  to  Figure  A.l  ( 1 )) 

2.  Next  procedure  is  that  the  f coordinate  system  is  rotating  to  the  g coordinate 
system  with  angle  0 with  respect  to  x-axis  (f, ).  (refer  to  Figure  A.l  (2)) 

3.  Final  procedure  is  that  the  g coordinate  system  is  rotating  to  the  b coordinate 

system  with  angle  y/  with  respect  to  z-axis  (g3).  (refer  to  Figure  A.  1 (3))  The  b 
coordinate  system  is  same  as  {P}  coordinate  system  in  Chapter  2. 


Figure  A.l  The  z-x-z  Euler  Angle  Expression  in  the  Coordinate  Systems 
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Figure  A.  1 can  be  regenerated  with  the  individual  coordinate  systems  as  shown  in 
figure  A.2. 


Figure  A. 2 The  z-x-z  Euler  Angle  Expression  in  Individual  Coordinate  Systems 


The  corresponding  expressions  for  the  rotational  matrices  for  the  individual 
coordinate  systems  are: 
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where 

e the  inertial  coordinate 

f the  first  intermediate  rotational  coordinate  of  first  Euler  angle  (j) 

g the  second  intermediate  rotational  coordinate  of  second  Euler  angle  0 

b the  third  intermediate  rotational  coordinate  of  third  Euler  angle  If/ 

~ Re  the  rotational  matrix  from  the  e coordinate  system  to  the  b coordinate 


system 


161 


CR^  the  rotational  matrix  from  the  b coordinate  system  to  the  e coordinate 


system 


For  finding  [U]  matrix,  the  angular  velocity  of  b coordinate  systems  relative  to  the 
e coordinate  systems  can  be  expressed  as  simply  adding  the  three  intermediate  angular 
velocities  expressed  in  the  same  coordinate  system. 
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where  5 to-  is  the  rate  of  change  in  the  b coordinate  system  with  respect  to  the  e 
coordinate  system.  The  scof  is  the  rate  of  change  in  the  f coordinate  system  with  respect 
to  the  e coordinate  system.  The  f to-  is  the  rate  of  change  in  the  g coordinate  system 


with  respect  to  the  f coordinate  system. 
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APPENDIX  B 

SUPPLEMENTAL  DERIVATION  DETAILS  OF  THE  EQUATIONS  OF  MOTION 
For  equation  2.55 
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The  following  are  the  expressions  for  each  of  the  terms  on  the  right-hand  side  of 
the  equation.  First,  for  the  platform  acceleration  dependent  terms,  the  forces  and 

moments  are  grouped  as 
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The  centripetal  acceleration-like  terms  are  grouped  in 
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The  gravitational  dependent  terms  for  the  platform  are  grouped  in 
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One  of  the  sources  of  nonlinearities  are  joint  frictions  which  are  denoted  by  the 
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The  friction  terms  F-  are  modeled  as  modified  friction. 
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The  left  hand  side  of  the  equation 
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The  derivation  of  the  Pscudo-Jacobian  matrix:  In  equation  (2.40),  the  F terms 
are  shown  on  F£  (equation  (2.42))  and  F/(  (equation  (2.43).  Thus,  the  Pseudo-Jacobian 
matrix  is  calculated  from  the  following  equations. 
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Therefore,  the  Pseudo-Jacobian  matrix  is  forming  with  (B.33)  and  (B.35)  as  shown 


on  equation  (B.21)  - (B.24) 


APPENDIX  C 

FIGURES  FOR  THE  COMPARISON  BETWEEN  FORCE-BASED  MODEL  AND 
ENERGY-BASED  MODEL  WITH  USING  THE  NOISY  DATA 
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mp  Value  (estimated  value)  - Force  without  Friction  - with  Noise  (0,  1 , 5 percent) 
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Figure  C.2  mp  Value  (0,  1 and  5 % Noise,  Force-Based  Model) 
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Figure  C.3  mb  Value  (0  and  1 % Noise,  Energy-Based  Model) 
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mp  Value  (estimated  value)  - Energy  without  Friction  - with  Noise  (0,  1 percent) 


Figure  C.4  mp  Value  (0  and  1 % Noise,  Energy-Based  Model) 
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Figure  C.5  me  Value  (0  % Noise,  Force-Based  vs.  Energy-Based  Model) 
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m Value  (estimated  value)  - Force  vs  Energy  - with  Noise  (0  percent) 
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Figure  C.6  mb  Value  (0  % Noise,  Force-Based  vs.  Energy-Based  Modeling) 
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Figure  C.7  mp  Value  (0  % Noise,  Force-  Based  vs.  Energy-Based  Modeling) 
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me  Value  (estimated  value)  - Force  vs  Energy  - with  Noise  (1  percent) 


FAXYWOWN  1 (Force) 

— EAXYWOWN  1 (Energy) 
True  Value 


i 

1 .5 1 1 1 1 __! L L 


o 5 10  15  20  25  30  35  40 

Iteration 


Figure  C.8  me  Value  (1  % Noise,  Force-  Based  vs.  Energy-Based  Modeling) 
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Figure  C.9  nib  Value  (1  % Noise,  Force-  Based  vs.  Energy-Based  Modeling) 
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mp  Value  (estimated  value)  - Force  vs  Energy  - with  Noise  (1  percent) 
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Figure  C.10  mp  Value  (1  % Noise,  Force-  Based  vs.  Energy-Based  Modeling) 
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